{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Case study.\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unrolling\n",
    "\n",
    "Let's simulate a kitten unrolling toilet paper.  As reference material, see [this video](http://modsimpy.com/kitten).\n",
    "\n",
    "The interactions of the kitten and the paper roll are complex.  To keep things simple, let's assume that the kitten pulls down on the free end of the roll with constant force.  Also, we will neglect the friction between the roll and the axle.  \n",
    "\n",
    "![](diagrams/kitten.png)\n",
    "\n",
    "This figure shows the paper roll with $r$, $F$, and $\\tau$.  As a vector quantity, the direction of $\\tau$ is into the page, but we only care about its magnitude for now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by loading the units we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "newton"
      ],
      "text/latex": [
       "$\\mathrm{newton}$"
      ],
      "text/plain": [
       "<Unit('newton')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radian = UNITS.radian\n",
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "N = UNITS.newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a few more parameters in the `Params` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Rmin</th>\n",
       "      <td>0.02 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rmax</th>\n",
       "      <td>0.055 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mcore</th>\n",
       "      <td>0.015 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mroll</th>\n",
       "      <td>0.215 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L</th>\n",
       "      <td>47 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tension</th>\n",
       "      <td>0.0002 newton</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>120 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Rmin           0.02 meter\n",
       "Rmax          0.055 meter\n",
       "Mcore      0.015 kilogram\n",
       "Mroll      0.215 kilogram\n",
       "L                47 meter\n",
       "tension     0.0002 newton\n",
       "t_end          120 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(Rmin = 0.02 * m,\n",
    "                Rmax = 0.055 * m,\n",
    "                Mcore = 15e-3 * kg,\n",
    "                Mroll = 215e-3 * kg,\n",
    "                L = 47 * m,\n",
    "                tension = 2e-4 * N,\n",
    "                t_end = 120 * s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` computes `rho_h`, which we'll need to compute moment of inertia, and `k`, which we'll use to compute `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params with Rmin, Rmax, Mcore, Mroll,\n",
    "                              L, tension, and t_end\n",
    "    \n",
    "    returns: System with init, k, rho_h, Rmin, Rmax,\n",
    "                         Mcore, Mroll, ts\n",
    "    \"\"\"\n",
    "    L, Rmax, Rmin = params.L, params.Rmax, params.Rmin\n",
    "    Mroll = params.Mroll\n",
    "    \n",
    "    init = State(theta = 0 * radian,\n",
    "                 omega = 0 * radian/s,\n",
    "                 y = L)\n",
    "    \n",
    "    area = pi * (Rmax**2 - Rmin**2)\n",
    "    rho_h = Mroll / area\n",
    "    k = (Rmax**2 - Rmin**2) / 2 / L / radian    \n",
    "    \n",
    "    return System(params, init=init, area=area, rho_h=rho_h, k=k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `make_system`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Rmin</th>\n",
       "      <td>0.02 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rmax</th>\n",
       "      <td>0.055 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mcore</th>\n",
       "      <td>0.015 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mroll</th>\n",
       "      <td>0.215 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L</th>\n",
       "      <td>47 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tension</th>\n",
       "      <td>0.0002 newton</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>120 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>theta               0 radian\n",
       "omega    0.0 radi...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>0.008246680715673206 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho_h</th>\n",
       "      <td>26.07109543981524 kilogram / meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>2.7925531914893616e-05 meter / radian</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Rmin                                              0.02 meter\n",
       "Rmax                                             0.055 meter\n",
       "Mcore                                         0.015 kilogram\n",
       "Mroll                                         0.215 kilogram\n",
       "L                                                   47 meter\n",
       "tension                                        0.0002 newton\n",
       "t_end                                             120 second\n",
       "init       theta               0 radian\n",
       "omega    0.0 radi...\n",
       "area                         0.008246680715673206 meter ** 2\n",
       "rho_h                26.07109543981524 kilogram / meter ** 2\n",
       "k                      2.7925531914893616e-05 meter / radian\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>theta</th>\n",
       "      <td>0 radian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>omega</th>\n",
       "      <td>0.0 radian / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>y</th>\n",
       "      <td>47 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "theta               0 radian\n",
       "omega    0.0 radian / second\n",
       "y                   47 meter\n",
       "dtype: object"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system.init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we compute `I` as a function of `r`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moment_of_inertia(r, system):\n",
    "    \"\"\"Moment of inertia for a roll of toilet paper.\n",
    "    \n",
    "    r: current radius of roll in meters\n",
    "    system: System object with Mcore, rho, Rmin, Rmax\n",
    "    \n",
    "    returns: moment of inertia in kg m**2\n",
    "    \"\"\"\n",
    "    Mcore, Rmin, rho_h = system.Mcore, system.Rmin, system.rho_h\n",
    "    \n",
    "    Icore = Mcore * Rmin**2   \n",
    "    Iroll = pi * rho_h / 2 * (r**4 - Rmin**4)\n",
    "    return Icore + Iroll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When `r` is `Rmin`, `I` is small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "6×10<sup>-6</sup> kilogram meter<sup>2</sup>"
      ],
      "text/latex": [
       "$6\\times 10^{-6}\\ \\mathrm{kilogram} \\cdot \\mathrm{meter}^{2}$"
      ],
      "text/plain": [
       "6e-06 <Unit('kilogram * meter ** 2')>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moment_of_inertia(system.Rmin, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As `r` increases, so does `I`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.00037418750000000006 kilogram meter<sup>2</sup>"
      ],
      "text/latex": [
       "$0.00037418750000000006\\ \\mathrm{kilogram} \\cdot \\mathrm{meter}^{2}$"
      ],
      "text/plain": [
       "0.00037418750000000006 <Unit('kilogram * meter ** 2')>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moment_of_inertia(system.Rmax, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "Write a slope function we can use to simulate this system.  Here are some suggestions and hints:\n",
    "\n",
    "* `r` is no longer part of the `State` object.  Instead, we compute `r` at each time step, based on the current value of `y`, using\n",
    "\n",
    "$y = \\frac{1}{2k} (r^2 - R_{min}^2)$\n",
    "\n",
    "* Angular velocity, `omega`, is no longer constant.  Instead, we compute torque, `tau`, and angular acceleration, `alpha`, at each time step.\n",
    "\n",
    "* I changed the definition of `theta` so positive values correspond to clockwise rotation, so `dydt = -r * omega`; that is, positive values of `omega` yield decreasing values of `y`, the amount of paper still on the roll.\n",
    "\n",
    "* Your slope function should return `omega`, `alpha`, and `dydt`, which are the derivatives of `theta`, `omega`, and `y`, respectively.\n",
    "\n",
    "* Because `r` changes over time, we have to compute moment of inertia, `I`, at each time step.\n",
    "\n",
    "That last point might be more of a problem than I have made it seem.  In the same way that $F = m a$ only applies when $m$ is constant, $\\tau = I \\alpha$ only applies when $I$ is constant.  When $I$ varies, we usually have to use a more general version of Newton's law.  However, I believe that in this example, mass and moment of inertia vary together in a way that makes the simple approach work out.  Not all of my collegues are convinced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object with theta, omega, y\n",
    "    t: time\n",
    "    system: System object with Rmin, k, Mcore, rho_h, tension\n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    theta, omega, y = state\n",
    "    k, Rmin, tension = system.k, system.Rmin, system.tension\n",
    "    \n",
    "    r = sqrt(2*k*y + Rmin**2)\n",
    "    I = moment_of_inertia(r, system)\n",
    "    tau = r * tension\n",
    "    alpha = tau / I\n",
    "    dydt = -r * omega\n",
    "    \n",
    "    return omega, alpha, dydt      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test `slope_func` with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0 <Unit('radian / second')>,\n",
       " 0.02939702689159846 <Unit('newton / kilogram / meter / radian ** 0.5')>,\n",
       " -0.0 <Unit('meter * radian ** 0.5 / second')>)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "slope_func(system.init, 0*s, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "results, details = run_ode_solver(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta</th>\n",
       "      <th>omega</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>115.2</th>\n",
       "      <td>205.83885219911747 radian</td>\n",
       "      <td>3.7827035471318 radian / second</td>\n",
       "      <td>36.27041961153422 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>116.4</th>\n",
       "      <td>210.40759585863975 radian</td>\n",
       "      <td>3.8320604323461227 radian / second</td>\n",
       "      <td>36.04569112433749 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117.6</th>\n",
       "      <td>215.0357989373311 radian</td>\n",
       "      <td>3.88180711854165 radian / second</td>\n",
       "      <td>35.818632262441746 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118.8</th>\n",
       "      <td>219.72393467957122 radian</td>\n",
       "      <td>3.9319528720968493 radian / second</td>\n",
       "      <td>35.589242955595836 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>120.0</th>\n",
       "      <td>224.47248760394635 radian</td>\n",
       "      <td>3.982507222010697 radian / second</td>\n",
       "      <td>35.357523131304944 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           theta                               omega  \\\n",
       "115.2  205.83885219911747 radian     3.7827035471318 radian / second   \n",
       "116.4  210.40759585863975 radian  3.8320604323461227 radian / second   \n",
       "117.6   215.0357989373311 radian    3.88180711854165 radian / second   \n",
       "118.8  219.72393467957122 radian  3.9319528720968493 radian / second   \n",
       "120.0  224.47248760394635 radian   3.982507222010697 radian / second   \n",
       "\n",
       "                              y  \n",
       "115.2   36.27041961153422 meter  \n",
       "116.4   36.04569112433749 meter  \n",
       "117.6  35.818632262441746 meter  \n",
       "118.8  35.589242955595836 meter  \n",
       "120.0  35.357523131304944 meter  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the results to see if they seem plausible:\n",
    "\n",
    "* The final value of `theta` should be about 220 radians.\n",
    "\n",
    "* The final value of `omega` should be near 4 radians/second, which is less one revolution per second, so that seems plausible.\n",
    "\n",
    "* The final value of `y` should be about 35 meters of paper left on the roll, which means the kitten pulls off 12 meters in two minutes.  That doesn't seem impossible, although it is based on a level of consistency and focus that is unlikely in a kitten.\n",
    "\n",
    "* Angular velocity, `omega`, should increase almost linearly at first, as constant force yields almost constant torque.  Then, as the radius decreases, the lever arm decreases, yielding lower torque, but moment of inertia decreases even more, yielding higher angular acceleration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot `theta`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xUVd7H8U86CSH00CGIcOgdRMSy2MsqdkUFrLgPru667tpWEctiX3XVB9til7Wj+Ii9gNJ7ywHpBEINpLeZef64QxyyJA6QzJ2ZfN+vV17J3Dsz93c0zDfn3nPPifH5fIiIiISbWLcLEBERORgFlIiIhCUFlIiIhCUFlIiIhKV4twuobcaYJGAQsA3wuFyOiIgcKA5oBcyz1pYE7oj6gMIJpxluFyEiItU6HpgZuKEuBNQ2gLfeeouWLVu6XYuIiATIzs7miiuuAP9ndaC6EFAegJYtW9K2bVu3axERkYP7r0swGiQhIiJhSQElIiJhSQElIiJhSQElIiJhqS4MkqhWbm4uO3bsoKyszO1Swl5CQgLp6emkpaW5XYqI1AF1OqByc3PZvn07bdq0ITk5mZiYGLdLCls+n4+ioiKysrIAFFIiUuvq9Cm+HTt20KZNG1JSUhROvyEmJoaUlBTatGnDjh073C5HRMJAYXEZm7fn1dr71+mAKisrIzk52e0yIkpycrJOh4oIhcVl/OXpH7np8e/IySuulWPU6YAC1HM6RPrvJSI+n4/n31/Klh35tEtPJS0lsVaOU+cDKpIUFxeza9euI36fzZs310A1IlJXfTlnIz8s2kK9xDhuHzWIuLjaiRIFVAS58sorWbx4MR9++CHnnXfeYb3Hd999x80331zDlYlIXbF+6z5e+GgZAOMu6kO7Fg1q7VgKqAiyZ8+eI36PnJwcvF5vDVQjInVNYXEZD782j7JyL6cd04GTBrSr1eMpoCLEuHHj2Lp1K7feeiu7d++mpKSE8ePHM3ToUIYNG8Z7771X8dzs7GzGjRvHMcccwymnnMKrr74KwNKlSxk/fjyrV69m4MCBAFhrueaaaxg2bBh9+vRh9OjRbN261Y0mikgY8/l8/OvdxWzdVUBGqzRuOL9XrR+zTt8HdTATXp7N/FXbQ3Ksgd1aMP66IUE997nnnmP48OHcdddd5Obmsn79ekaOHMl9993H1KlTufvuuznjjDNISUnhxhtvZODAgfz4449s27aNsWPH0qhRI0aMGMGECRN47bXXmDp1KgA333wzF198MS+//DL5+fmMGzeOl156ifHjx9dm00Ukwnz203pmLtlKclI8d4weRFJCXK0fUz2oCNWqVStGjRpFTEwMZ511FuXl5WRnZ7N8+XI2bdrEHXfcQVJSEhkZGVx99dVMmTLloO/z0ksvMWbMGMrKysjOzqZx48a6z0lEDrB6Uw6vfLIcgJsv7Uub5qkhOa56UJUE26NxW8OGDSt+Tkx0hniWl5eTlZVFUVERQ4b82g6v10ujRo0O+j7Lly9n7Nix5OXl0blzZ4qKimjSpEntFi8iESOvsJRHXp9HucfHOcM6MqxPm5AdWwEVZdLT02natCkzZ/66cvKePXsoLv7vG+m2b9/Obbfdxptvvkn//v0BePDBB3UNSkQA8Hp9PPn2QnbkFNGlfSOu+X3PkB5fp/giSEJCAnl51U8r0rt3b1JTU3n++ecpLS1lz549/M///A/PPPMM4PS2CgoK8Hq95Ofn4/P5qFevHgCzZs1i6tSpmilCRAB475vVzF+1nQYpCdx+1SAS4kMbGepBRZALLriACRMmkJqaStOmTQ/6nMTERF588UX+8Y9/cPzxxxMTE8Mpp5zCXXfdBcCgQYOIj49nwIAB/PDDD9x8881ce+21lJeX07FjRy6//HKmTZuGz+fTrBEiddgiu4O3vsgkJgb+csUA0pukhLyGGJ/PF/KDhpIxJgNY/80339C2bdsD9q1atYpu3bq5Ulck0383kei2M6eIP/3ze3ILSrn8NMPI07vW2rG2bNnCySefDNDRWrshcJ9O8YmISIWycg+PvD6P3IJS+nVpzqWnGtdqUUCJiEiFl6Yux27KoXnjZP5yxQDiYt071a+AEhERAL6dv4nPf95AfFwsd44eRMPUJFfrUUCJiAjrsvbx3HtLALjxgt50btfY5YoUUET7IJGapv9eItEnr7CUia/NpbTcy6mD23P6kA5ulwTU8YBKSEigqKjI7TIiSlFREQkJCW6XISI1xOP18fhbC8jeXUintg258YLebpdUoU4HVHp6OllZWRQWFqpn8Bt8Ph+FhYVkZWWRnp7udjkiUkPe+SKThZk7aJCSyF2jB5MYgklggxXSG3WNMacCDwOdgR3AY9baF4wxicCzwEWAB3jSWjsx4HWXAP8AWgE/AGOstUc8o2laWhoAW7du1ewJQUhISKBFixYV/91EJLLNWraN/3y9mtgY+NtV7tyMW52QBZQxph3wATAamAoMAL4wxmwATgIM0AloCEw3xmRZa183xnQHXgHOBOYDjwBTgOE1UVdaWpo+cEWkztm8PY9/vrMQgNFnd6dvl/A7MxLKU3wZwNvW2o+stV5r7Tzge+A4nNB6yFqb47+T+HFgrP91VwKfWmtnWmuLgTuB44wxnUNYu4hI1CgoKuOhyXMpKinnuN6tOf+ko90u6aBCFlDW2hnW2hv3PzbGNAGOBxbhnLpbGfD0TGD/co3dA/dZawuBzQH7RUQkSPtnKM/amU9GqzRuuaxf2M676cogCWNMQ+ATYA6wwL+5MOAphcD+k6GplfZV3i8iIkF650vL3JXZpCYncNeYwSQnhe+c4SEPKGNMF2A2sB1nUMT+9SOSA56WAuT7fy6otK/yfhERCcKsZduY8pUlNgb+euVAWjWr73ZJ1QppQBljTsDpNX0MXGStLbbW5gDZOIMk9uvKr6f1VgbuM8akAO058JSgiIhUY+O2XP75jnPCatRZ3enfNfwGRVQWylF8nYBpwN3W2n9V2v0GMN4YsxTnlN5twNP+fW8DM40xJwGzgInAImvt6pAULiIS4fIKS3lw8hyKSjyc0K8NF/wuPAdFVBbKk4/jgAbARGPMxIDtzwH3Ak8AK3B6dS8CkwCstcuMMdf4H7fB6YFdHMK6RUQilsfj5dE35lfMFPHHS/qG7aCIykIWUNbaW4Fbq3nKOP/XwV77Ac49VCIicggmT1vJ4tU7aZSaxF1jBlMvMXwHRVRWp6c6EhGJZl/N2cjUH9cSHxfDHaMHkd44sgY/K6BERKLQyvW7ef6D/ctn9KHHUU1drujQKaBERKLMjpxCJr46j3KPj3OGdQyb5TMOlQJKRCSKFJeU89C/57I3v4S+nZtz3bk93S7psCmgRESihNfr48l3FrJu6z5aNavP30YNJC4ucj/mI7dyERE5wJvTVzFr2Tbq14vnnmuOoUFKotslHREFlIhIFPh+wWbe+2YNsbEx/G3UINq1aOB2SUdMASUiEuEyN+zhmXcXA3D9eT3pb8J/GqNgKKBERCLY9j2FPDR5LmXlXs48NoOzj+vodkk1RgElIhKhCovLeOCV2c6IvS7NueH8XhEzjVEwFFAiIhHI4/HyyBvz2ZidR7sWqdw+ahDxETxi72CiqzUiInXEy1OXszBzBw1SErn32iGkJie4XVKNU0CJiESYT2asZdpP64mPi+XuqwfTsml4Lzx4uBRQIiIRZO6KbF6euhyAWy7rF5Fz7AVLASUiEiHWbtnLY2/Ox+eDkad35aT+bd0uqVYpoEREIsDOnCLuf2UOxaUefjegLZed2sXtkmqdAkpEJMwVFJVx/yuz2ZNbTI+jmkbUqrhHQgElIhLGyj1eHn59Hhu25dI2PZW/Xz2YhPg4t8sKCQWUiEiY8vl8PP/+kool28dfN4TUCJ8A9lAooEREwtR/vl7NV3M3kZgQxz3XHhO1w8mrooASEQlDX8/dxFvTM4mJgduuGECX9o3dLinkFFAiImFmod3Bs+85s5PfMKIXx/Zq5XJF7lBAiYiEkXVZ+3j4tXl4vD4uOOlozhl2lNsluUYBJSISJrbvKeS+l2ZRVFLO8X3bMPrs7m6X5CoFlIhIGMgtKGX8i7PIySuhV6dm/PnyfsTGRv+9TtVRQImIuKy4tJwHXplN1s58MlqlcXcdutepOgooEREXeTxeHn9zAZkbc2jeOJn7rh9C/ShcOuNwKKBERFzi8/l47v0lzFmRTWpyAvddN4SmDZPdLitsKKBERFzy5vTMihtxx183hPYt09wuKawooEREXDBt5jre/Xo1sbEx3D5qIF0zmrhdUthRQImIhNiPi7bw4sfLAPjjxX0Y3L2lyxWFJwWUiEgILbQ7+Oc7C/H5YNRZ3ThlcAe3SwpbCigRkRCxG/cw8dW5lHt8jDixExcN7+x2SWFNASUiEgKbsnOZ8PJsiks9DB/YjqvP6VEnFh08EgooEZFatn1PIfe+OIu8wjIGdW/BHy/pW+dniQhGvBsHNcYMBqZZa9P9j5OAPKA04Gk/W2tP8++/BPgH0Ar4ARhjrd0R2qpFRA5dTm4x90z6md37nOXabx81iPg49Q2CEdKAMsbEANcCj1fa1QvYY639r6EsxpjuwCvAmcB84BFgCjC8dqsVETky+YWl3PviLLbtLqBT24bcc80xJCVoCqNghTrGJwB/AB6stH0AsLiK11wJfGqtnWmtLQbuBI4zxujqooiEraKScu57eTYbtuXSpnkqE64/VlMYHaJQB9Qka+0AnJ5QoP5AujFmqTFmuzHmPWNMG/++7sDK/U+01hYCm3F6XSIiYae0zMNDk+dg/fPr3T/2WBqmJrldVsQJaUBZa7dWsasA+Ak4GTBAEfCRf18qUFjp+YVASm3UKCJyJMo9Xh59Yz5L1uyiUYMkHhw7lPTG+rg6HEFdgzLGNMUJj4FAOuABsoF5wFfW2oIjKcJae2ul490K7DTGtMMJr8qzJ6YA+UdyTBGRmubx+njqnUUVk78+MHYorZunul1WxKq2B2WMOdoYMxnIAp4A+gCJOL2aIcBLwC5jzEtHck3IGHO/MaZbwKZE//dinNN7JuC5KUB7Ak77iYi4zefz8fz7S/hh0RaSk+KYcMOxZLTS5K9HosoelDHmXmAU8BowwFq7oorn9QCuAL4yxrxqrb3vMOroDQw0xoz0P34a+Mxau9MY8zYw0xhzEjALmAgsstauPozjiIjUOJ/Px0tTl/PlnI0kJsRxz7VD6NK+sdtlRbzqelAFQHdr7QNVhROAtXaFtfYuoBuHf9rtWiAH+AXYgHM/1FX+918GXANMAnYBPYCLD/M4IiI1yufz8cbnq/h0xjri42K5e8xgenVq5nZZUSHG5/O5XUOtMsZkAOu/+eYb2rZt63Y5IhJl3vnS8vYXmcTGxnDn6EEM6dnK7ZIiypYtWzj55JMBOlprNwTuq+4U36hgD2Ctff2wqxMRiVDvf7vGCacYuG3kAIVTDatuFN8TlR43AbzAVpxTcO2AOGANoIASkTpl6o9ree2zlcTEwC2X9ef4fm1++0VySKoMKGtt8/0/G2PGARcAo6y1Wf5tzYDJwMLaLlJEJJx89tN6Xp66HIBxF/Vl+MB2LlcUnYK9UXc8cMv+cAKw1u4C7gBuro3CRETC0eezNjDpw6UA3Hh+L04fogUHa8uhzCRxsJOrnXDuVRIRiXpfztnI8+8vAeD683py9rCjXK4ougU7m/lk4DVjzP3AIiAGOAa4G3iqlmoTEQkbX8/dyLPvOXNaX3tuD849oZPLFUW/YAPqTpz58cYDLfzbtgGPWmsrL50hIhJVvp67kWfeXYzPB2PO7s6IE492u6Q6IaiAstZ6gfuA+/yDI/ZfgxIRiWqB4TTqrG5cOFwr/YRK0AsWGmP64Sx9Eed/HAMk4UyDNLZ2yhMRcU/lcLr45C5ul1SnBDub+d3AAzhTGdUH9gEN/bv/r3ZKExFxzxezN/Lc+wonNwU7im8s8FdrbRrOtafeQBtgNs6SGyIiUePzWRt49r1frzkpnNwRbEC1BD7w/7wYONZamw38Df+kriIi0eCzmesqhpJfe24PXXNyUbABtRNo6v95Nc66UOCsE9W6posSEXHD1B/XMumjZQBcP6KnRuu5LNiAmgq8aIzpC3wHjDLGnAjcCmysreJERELl/W/XVExfdOP5vTj3eN3n5LZgR/HdBjwJ9ATeAi4EvsEZNDGymteJiIQ1n8/HlK9W8/YXmcTEOHPrafqi8BBsQF0I3G2t3e1/PMYYcxNQbK0tr53SRERq1/7FBt/7Zg2x/lnJNfFr+Ag2oJ4B5gL7Awpr7eGunisi4jqv18fLnyzn0xnriIuN4S8jB2jJjDAT7DWoOcD5tVmIiEioeLw+nn1vccUy7XeOHqRwCkPB9qC8wD+MMX8H1uPMy1fBWju4pgsTEakN5R4v/3xnIT8uyiIxIY6/Xz2Yfibd7bLkIIINqDn+LxGRiFVS5uGR1+cxb+V2kpPiuffaY+jZqZnbZUkVqgwoY0yctdYDYK2dEMybBb5GRCScFBaX8dDkuSz9ZRcNUhKYcMOxdG7X2O2ypBrVXYOaZ4y5xD8pbLWMMfHGmCuB+TVXmohIzcgtKOWeF35m6S+7aJKWxMRxwxROEaC6U3znAf8CnjbGTAWmAyuAXTgLFjbHmVHiROAinFF+I2q1WhGRQ7R7XxH3vDCLzdvzSG+SwkM3DqVl0/pulyVBqDKgrLWbgRHGmP7AH4FJQDrgC3jaduBz4GxrrXpPIhJWtu7M554XfmZHThHtWzbg/huOpWnDZLfLkiD95iAJa+1C4GoAY0x7nBV1vUC2tTardssTETk867L2Mf7FWezNL8G0b8z464fQICXR7bLkEAS9YCGAtXYTsKmWahERqRHLftnFg5PnUFhcTt/Ozbnr6sEkJx3Sx52EAf0fE5GoMmvZVh57cwFl5V6G9WnNrSP7kxAf53ZZchgUUCISNb6YvYHn31+C1wdnDc3ghvN7Exf7mwORJUwpoEQk4gXOSA5w+WmGy08zxMQonCLZIQWUMSYV6AysBBKttXm1UpWISJA8Xh+TPlzK9FkbiI2BGy/sw5nHZrhdltSAoALKGJMIPAXc4N/UBXjEGJMMXGGt3VdL9YmIVKmkzMPjb85n9vJsEuJj+euVAzm2Vyu3y5IaEuxs5g8AQ4HjgWL/tseADJyFDEVEQmpffgl//9+fmL08m/rJCTwwdqjCKcoEG1CXADdZa2fhv1HXWjsXuB74fS3VJiJyUNm7C7j92RlkbsyheeNkHr1pGD2Oaup2WVLDgr0GlQ5kH2R7LpBSc+WIiFRvzeYc7n9lDnvzSujYOo3x1w3R7BBRKtge1I/ALQGPff7rUvcAM2u8KhGRg5izfBt3Pv8Te/NK6Nu5OQ+PG6ZwimLB9qBuBr4wxpwK1ANexRnN5wFOq53SRER+NW3mOl76eBleH5w8qB3jLupLQnywf2NLJAoqoKy1a4wx3YCRQHf/694C3rTWFh7qQY0xg4Fp1tp0/+NE4FmcWdE9wJPW2okBz78E+AfQCvgBGGOt3XGoxxWRyOPx+pj86Qqm/rgWgJGnd+WyU7voHqc6IOj7oKy1JcDkIzmYf22pa4HHK+2aABigE9AQmG6MybLWvm6M6Q68ApyJs97UI8AUYPiR1CIi4a+4pJzH31rAnBXZxMfF8MdL+jJ8YHu3y5IQqW5F3XkcuLRGlay1g4M83gTgbOBB4O8B20fj9IpygBxjzOPAWOB14ErgU2vtTH9dd/qf09lauybI44pIhNm9r4gH/j2HtVv2kZqcwF1jBtPraC3PXpdU14OaVgvHm2StvdcYc9L+DcaYRjin7lYGPC8T6OX/uTsBK/VaawuNMZv9+xVQIlFo7Za9PPjvOezaV0yrpvW597pjaJvewO2yJMSqW7BwQk0fzFq79SCbU/3fA69lFfLr8PXUSvsq7xeRKDJ7+TYef2sBJaUeumU04e6rB9MwNcntssQFwU519O8qdvmAUiAL+MBau+owaijwfw8cK5oC5AfsrzyONHC/iEQBn8/HR9//wqufrcTng+ED23HTxX20VEYdFuwYzTyc60TdgL3+r844K+22AIYA840xZxxqAf7rTtk4gyT268qvp/xWBu4zxqQA7TnwlKCIRLCycg9PTVnE5GlOOI06qxt/uqyfwqmOC3YU31HAw9bauwM3GmPuAfpba88xxtyAM/hh+mHU8QYw3hizFOeU3m3A0/59bwMz/detZgETgUXW2tWHcRwRCTM5ecVMfHUeqzbsISkxjj9f3p/jerd2uywJA8H2oIbj3Jxb2RTgdP/P03F6WIfjXmA5sAKYB3wATAKw1i4DrvE/3gX0AC4+zOOISBhZl7WPvzz9I6s27KFZo2QeGTdM4SQVgu1BbcYJosqj5s7g1zn6OgA5wbyZtfZ7oFHA42JgnP/rYM//ACe0RCRKzFicxVNTFlFa5sF0aMzdYwbTOK2e22VJGAk2oO4F3vSfZpuH0/MagDOT+dX+m2nfAt6pjSJFJHp4vT7e+iKTd792ztIPH9iOcRf1ITFB15vkQEGd4rPWvgucBJTg3Dh7Mc5Q76HW2rdxrhs9CtxeO2WKSDQoKCrjoclzeffr1cTGwHXn9eRPl/VTOMlBHcpURz8DP1exby4wt6aKEpHos3l7Hg9NnkvWznxSkxP461UD6W/S3S5Lwliw90ElAzfinNZLAA6YpdFae0nNlyYi0WLO8m088fZCikrKyWiVxl1jBtOqWX23y5IwF2wP6kXgApyRerm1V46IRBOP18c7X2TyH//1pmF9WnPLpf2olxT0yRupw4L9LTkDGGmtnVqbxYhI9MgrLOXxtxawMHMHsTFw1VndufB3R2uZDAlasAFVBujGWBEJyrqsfUx8bS7ZuwtpkJLI7VcNpE+X5m6XJREm2Bt1/wk8bIzRb5iIVOvruRv56zM/kr27kKPbNuSpP5+ocJLDEmwP6hKgN5BtjMnDmSC2wv6VcUWk7iot8/DCR8v4cs5GAE4f0oEbRvTSEHI5bMEG1LNVbG9KpbASkbpn264CHn5tHuu27iMxPpY/XNibUwZ3cLssiXBBBZS19rXAx8aY03Dmxxvhf4+qAkxEotzPS7fy9H8WUVhcTqum9bl91EA6tW302y8U+Q1Bj/U0xmTgLK8xBmiLsx7TiyicROqksnIvr362gk9+XAfAsb1accul/aifnOByZRItqg0oY0wScBFOb+lEwAt8D7QBTrDWLqntAkUk/GTvLuDRN+azZvNe4mJjuPr3PTj3+KM0hFxqVJUBZYx5HrgcSAS+Aq4FPrHW5hhjynCGnotIHTNr2VaenrKIguJy0hsn87erBmI6NHG7LIlC1fWgbsS59+kh4P+stbtDU5KIhKPSMg///nQFn/20HoBjerTkT5f1IzUl0eXKJFpVF1C/A67AWdn238aYn3DWZPooFIWJSPjYvD2PR9+Yz4ZtucTHxTDmHJ3Sk9pX5Y261tofrLU3AC2By4A9wGPARv/rLjDGpIWkShFxhc/n48s5G/nzUz+wYVsurZrV57E/nsB5J3RSOEmt+81RfNbaUpye0wfGmEbApTg9qwnAHcaYKdba62q3TBEJtfzCUp59fwk/LdkKwEn92/KHC3uTUk+j9CQ0DmlKYWvtXuAF4AVjTHvgKmBkbRQmIu5ZsW43T7y9gJ05RSQnxfGHC/vwuwHt3C5L6pjDnvPeWrsJZwDFQzVXjoi4qdzj5e0vMvng2zV4fdClfSNuu2Kg1m4SV2hRFhEBIGtnPk+8tYA1m/cSGwOXnNKFy08zxMcFO6e0SM1SQInUcT6fj89nbeCVT1ZQWuYhvXEyt44cQI+jmrpdmtRxCiiROmz3viKeeXcxCzN3AHDSgLbceH5vTVckYUEBJVJHzViUxf9+uIS8wjIapCTwPxf1YVifNm6XJVJBASVSx+zLL2HSh0uZ6R8+3t+kc/OlfWnaMNnlykQOpIASqUPmLN/Gc+8vISevhHqJcVxzbk/OGNJBN91KWFJAidQBeYWlvPjxMr5fsAWA7h2b8KfL+mv4uIQ1BZRIlJu7Ipvn3l/MntwSEuNjueqs7vz++KOIi1WvScKbAkokSu3LL+HFj5fx46IsALplNOFPl/WjdfNUlysTCY4CSiTK+Hw+Zi7eygsfL2VffimJCXGMOqsb5wxTr0kiiwJKJIrs2lvE/36wlLkrswHofXQzbrq4r641SURSQIlEAa/XmQ3itc9WUlRSTkq9eMac04PTj+lArHpNEqEUUCIRbuO2XJ59bzGZG3MAZ6XbP1zYW/c1ScRTQIlEqJIyD//5yvLhd7/g8fpokpbEDSN6M7R3K93XJFFBASUSgRZkbmfSh0vJ3l1ITAycOTSD0Wd11xx6ElUUUCIRZPe+Il6aurxildsOLRsw7qK+dOvYxOXKRGpe2ASUMeYanNV6SwI2jwPeAZ4FLgI8wJPW2omhr1DEPeUeL5/OWMc7X2ZSVOIhKTGOkad15dwTjtJ6TRK1wiaggP7AE9baOwI3GmMmAgboBDQEphtjsqy1r7tQo0jILVu7i0kfLmVTdh4AQ3q25PoRvUhvnOJyZSK1K5wCagDw9EG2jwbGWGtzgBxjzOPAWEABJVFt194iJn+6gh8XOzNBtGpanxvO78XAbi1crkwkNMIioIwxcUBv4CpjzJNAIfAyzim/VsDKgKdnAr1CXqRIiJSVe/j4h7W8+/Vqiks9JMbHctHwzlw4vDOJCXFulycSMmERUEBzYD7wGnAB0A2YCiT69xcGPLcQ0LkNiTo+n485K7L59ycr2La7AIBje7Xi2nN70qKJfuWl7gmLgLLWZgMnBmxabIz5F3Cm/3HgHYcpQH6oahMJhY3bcnn5k+UsXr0TgHYtGnD9eT3pZ9JdrkzEPWERUMaYHsAl1trxAZsTgWIgG2eQRJZ/e1cOPOUnErH25pXw1heZfDl7A14fpCYncMUZXTnz2AziNDpP6riwCChgL/AXY8wW4BWgH3AzcBOwAhhvjFkKpAK3cfDBFCIRo6TMw6cz1vHu16spKiknNjaGs4dmMPL0rqTVT/ztNxCpA8IioKy1WcaYc4FHgX8Cu4AHrLXvG2OmAUdv2NYAABC3SURBVE/gBFUs8CIwybViRY6A1+vjh0VbeP3/VrFrbxEAA7u14Jrf96BdiwYuVycSXsIioACstd8CAw+yvRjnht1xIS9KpAYtXr2DVz9bydot+wDo2DqNq8/poetMIlUIm4ASiVa/bNnLa5+trBgA0bRhPa46sxsnDWinBQRFqqGAEqklWTvzeWt6JjP8N9rWrxfPRSd34ZxhHamXqH96Ir9F/0pEatiuvUVM+cry1dxNeL0+EuJjOfu4jlxyShcapGgAhEiwFFAiNSQnt5j3vl3D9FkbKCv3Ehsbw2nHdOCyUw3NG2vxQJFDpYASOUL78kv48LtfmPbTekrLPAAc37cNV5zRlTbNU12uTiRyKaBEDtO+/BI++t4JppJSJ5iG9GzJyNO70rF1Q5erE4l8CiiRQ5STV8xH36/l85/XU+wPpkHdW3D5aYbO7Rq7XJ1I9FBAiQRp974iPvzuF6bP2kBpuRdwbrK9/DRDl/YKJpGapoAS+Q1bd+Xz4Xe/8M28TZR7fAAc06Mll51qOLpdI5erE4leCiiRKqzdspcPv/uFmUuy8PogJgaO69OaS0/pomtMIiGggBIJ4PP5WLpmF+9/t6Zi5oe42BhOHdSOC4d31qg8kRBSQIkA5R4vMxZn8fH3a1m31Zkrr15iHKcPyeC8EzrpPiYRFyigpE7LLyxl+uyNfDZzHbv2FQPQKDWJc4Z15KzjOmrmBxEXKaCkTtqyI49PZ6zjm/mbK+5hatcilREnHs1J/duSmBDncoUiooCSOsPj9bEgczvTZqxjkf/6EkC/Ls0594RO9DfpxGp2cZGwoYCSqLcvv4Rv5m3i81kbyN5dCEBiQhy/G9CW3x9/FB1aprlboIgclAJKopLP5yNzQw7TZ29gxuIsyvw31qY3SeHsoR059Zj2ur4kEuYUUBJV8ovK+GHBZqbP3siGbbmAc//SwG4tOGtoBv27ttAigSIRQgElEc/n87Fi3W6+nLORn5ZsrZiGqGFqIqcMas8Zx2bQsml9l6sUkUOlgJKItTOniG8XbOKbeZvZtqugYnvvo5txxrEZDOnZioT4WBcrFJEjoYCSiFJcUs6s5dv4dv5mlqzZic+ZGo8mafU4ZXB7Th3cXr0lkSihgJKw5/F4WfLLLn5YuIWfl26tWOIiPi6WIT1bcvKg9vTr0py4OPWWRKKJAkrCks/nw27K4cdFWcxYnMXevJKKfV07NGb4wHYM69tGI/FEopgCSsKGz+dj7ZZ9zFzihNKOnKKKfa2b1efE/m05aUBbWjfThK0idYECSlzl9fpYszmHn5du46elW9m+p7BiX9OG9RjWpw0n9m/D0W0bEROj4eEidYkCSkKu3ONlxdrdzFq+jdnLt7HbP0krQOMGSQzt3Zrj+7ahW0YTTT0kUocpoCQk8ovKWJi5nTkrslmwajsFxeUV+5o1rMexvVtzXO/WdM1oohtpRQRQQEkt8fl8bNqex/yV25m3ajurNuzB6/VV7G/XIpUhPVsxpGcrOrfT6TsR+W8KKKkxeYWlLF2ziwWZ21lkd1SsrwTOqrS9OjVjcI+WDO7RQgMdROQ3KaDksJWUechcv4clv+xk0eqdrN2yt+LGWYBGDZLob9IZ1L0F/bqkUz85wb1iRSTiKKAkaCVlHuzGPSxfu5tla3eRuSGHco+3Yn98XAzdMprSzzRnQNcWZLRK0yAHETlsCiipUm5BKZkb9rBy/W5Wrt/Dms05lHt+7SLFxMBRbRrS++hm9O3SnB4dm1IvSb9SIlIz9GkigDP0e8O2XNZsyiFzYw524x6ydhYc8JyYGDiqdUN6dmpKj6Oa0uvoZprJQURqjQKqDvJ4vGzansfaLftYu2Uva7bsZX3WvoplKvZLjI+lc/vGdO/YhG4ZzleqAklEQkQBFeX25ZewMTuXDdty2bA1l/Vb97ExO69ihdlArZvVp3O7xnTNaIzp0JiOrRsSrwlYRcQlCqgo4PX62L2vmKydeWzZkc/m7c73Tdl57M0vOehrWjZNoVPbRnRq05Cj2zaic7tG6h2JSFiJiIAyxvQBJgG9gXXANdbaee5WFVolZR525hSyY08R2/cUsG13Idm7C9i2q4CtuwooLfMc9HXJSXG0b5lGh5ZpdGydRsfWDclolaYh3yIS9sI+oIwxicBU4CngBOBC4EtjTAdrba6rxdUAj8dLbkEpe/NL2JtXwp7cYudrXzE79xaxa18Ru/cWV9kT2q9RgyTaNE+lTfNU2rVIpW16A9q1aEB642TN0iAiESnsAwo4CUiw1j7lfzzFGHMTcCnwUigK8Hi8OLP0+PD6nMcer49yj5fych9l5R5Ky72UlnkoLi2npNRDcYmHwpJyikrKKCwup6CojPyiMgqKysgtKK34yi8qPeDm1qrExcbQvHEy6Y1TaNEkhZZN69OqaX1aNkuhdbNU9YhEJOpEQkB1B1ZV2pYJ9ArFwd+cvor/fLW6Vo/RICWRRg2SaNwgiSZp9WiSVo/GafVo1qgezRol07xRMo0a1NMkqiJSp0RCQKUChZW2FQIpoTh4fFws8XEx+HwQExNDTIwzY0JsbCxxsTEkxseSEB9HQkIsifGxJCXGUy8xjqTEOFKSEkipF09yUjypKQmkJidQPzmBBimJpNVPpEH9RBqkJGqknIjIQURCQBUAyZW2pQD5oTj4ZacaLjvVhOJQIiISIBL+dF8JVE6Irv7tIiISpSKhB/UdEGOM+TPwLM4ovt7AR65WJSIitSrse1DW2lLgTJxg2gPcDYyw1u50tTAREalVkdCDwlq7HBjmdh0iIhI6Yd+DEhGRukkBJSIiYUkBJSIiYSkirkEdoTiA7Oxst+sQEZFKAj6b4yrvqwsB1QrgiiuucLsOERGpWitgbeCGuhBQ84DjgW3AwdekEBERt8ThhNN/LaEU4wtmKm0REZEQ0yAJEREJSwooEREJSwooEREJSwooEREJSwooEREJSwooEREJSwooEREJSwooEREJS3VhJonDZozpA0zCWcF3HXCNtfa/7naOBMaYU4GHgc7ADuAxa+0LxphEnJWKL8KZaeNJa+1E9yo9PMaYRsBS4F5r7avR0C5jTCvgf4HfAcXAi9bae6KkbUOAZwAD7AQetta+HMltM8YMBqZZa9P9j6ttizHmEuAfOLMo/ACMsdbuCHnhQThI29KBp4GTgRjgc+AWa22Of3+NtE09qCr4f7mmAv8BGgEPAV8aY9JcLewwGGPaAR8AD+K05XJgojHmdGACzodEJ2AQMNoYM8qtWo/AJKBNwONoaNdUnCm6WgBDcNowkghvmzEmFqdtz1hrG+L8Pj7r/4Mw4tpmjIkxxlwHfAkkBuyqsi3GmO7AK8AYoCmwBpgSwrKDUk3bXgbKgY44f/Q2Bp7zv6bG2qaAqtpJQIK19ilrbZm1dgqwArjU3bIOSwbwtrX2I2ut198L/B44DhgNPGStzbHWbgAeB8a6VejhMMaMBtKAZQGbI7pdxphjgKOAm621xdba9Ti/k98R4W3D+TBLB2KMMTGAD+fDrpTIbNsE4A84fwAGqq4tVwKfWmtnWmuLgTuB44wxnUNUc7D+q23+PzC8wARrbYG1di/wEr+uel5jbVNAVa07sKrStkyglwu1HBFr7Qxr7Y37HxtjmuBMoLsIpwu+MuDpEdVGY0xHYDxwTcC2RkR4u4ABOIF7nzEmyxizFjgfKCLC22at3Y1z6us1oAxnktC7cHqLkdi2SdbaAcD8/RuC+B3sHrjPWlsIbCb82vpfbfP/kTvCWvtLwPNG4HyeQA22TdegqpYKFFbaVgikuFBLjTHGNAQ+AeYAC/ybA9sZMW00xsQBbwK3WWuzjTH7d6X6v0dku/z2/xHxA05PqiswHed6DURw2/x/gRcDI3FOPQ8FPgT2+p8SUW2z1m49yObf+h2MiM+XKtp2AGPMbTgBNdS/qcbapoCqWgGQXGlbCpDvQi01whjTBefc/0rgCn5tX2A7I6mN9wDWWvthpe0F/u+R2i6AEiDXWnuf//ESY8zLOKeNILLbdgFwnLX2r/7HPxhjXiE62rbfb/0ORvznizEmAfgX8HtguLU207+rxtqmU3xVW4lzgTNQVw7sskcMY8wJOL2mj4GL/Nc1coBsDmxnJLXxMuAiY8xeY8xenFMIz+MMaInkdoFzOijFP1hnv3gg0v+fAbQDkiptK8fpHUZ62wAI4t/WAZ8vxpgUoD0R0lZjTAPgK5zBH4OttYsDdtdY29SDqtp3OBdx/4xzvvxCnOHmH7la1WEwxnQCpgF3W2v/VWn3G8B4Y8xSnK75bTjDR8OetbZr4GNjzGLgKf8w83witF1+X+F8YD9hjPkLzj/4a3EuWK8jstv2Jc4o0htwLq73B64HrgM2EdltC1Tdv623gZnGmJOAWcBEYJG1drUbhR6GKTgdnOP915gC1Vjb1IOqgrW2FDgTJ5j2AHcDI6y1O6t9YXgaBzTA+VDID/h6BLgXWI4zQnEezjWBSe6VWmMiul3+0U8n4lx/2oZz/elRa+0HRH7bVuCc5huLc93pbeAOa+1UIrxtlVTZFmvtMpyBPZOAXUAP4GJ3yjw0xpjewFnAYGBHwOfJFqjZtmlFXRERCUvqQYmISFhSQImISFhSQImISFhSQImISFhSQImISFhSQImISFjSjboiNcQY8yq/TtdzMBNwZpH/DmhgrQ3JtDb+OQt/AkZVd7Okf4682cBV1lobitpEqqMelEjNuQVnButWOEtjgHMz4/5tjwM/+38uOMjra8vNwJLfupPfWusF7idyb4yVKKMbdUVqgTGmJ85yGR39awG5VUc9nOmDhltrlwf5mrXAtdba72uzNpHfolN8IiHkn5+s4hSfMcaHs6LsnTjz7c3HWfDtr8BVQC5wp7X2Df/rGwBP4Cwj7gO+xVlqu6plES4D9gaGkzHmHuAGoDnOmmd3WWs/D3jNRzi9we9roMkih02n+ETc9zDwJ5xl3dsDC3GCaRDOOkkvGGP2ry/0Ik6QnY4zV58P+MIYU9Ufm2fjzOMHgDHmfP+xrsSZXfsz4D1jTFrAa6YDp1TzniIhoYAScd9z1trv/EsWTMNZN+cu/0CFJ3HW1ulojDkKp0c00lo7z98rugrIAM6o4r0H4kxWul8GzlpTG/2nHu/Hmbi1LOA5K3Fm3z5gtniRUNNfSCLuC1w6uxDYYK3df3G42P89Cejg/9kGrB4MzmJwBifcKmuBM6P0fm/ijDRcZ4xZgLO68mRrbVHAc3b7v6cfYjtEapR6UCLuK6v02FvF8+L9z+0H9A346gJMruI1XiBm/wP/cjEDcHpcPwNjgKX+QR377f9c8ATdApFaoIASiRyrgASgvrX2F2vtLzhrRT2GE1IHk40zGAIAY8wFwFhr7ZfW2ltwel55OOv77Nc84LUirtEpPpEIYa21xphPgNeNMeNwVtx9CGdwRWYVL1sA9Al4HAc8ZozZjjNicAjQ0v/zfn1wlpYPPPUoEnLqQYlEltE4YfIxziqtDYFTrbV7q3j+Zzij/QCw1r4HjMfpda0GHgRustZ+G/CaE4Dp1lqd4hNX6UZdkShmjEkBNgBnWGsXBvH8WGAjzkjBGbVcnki11IMSiWLW2kKc3tK4IF9yHrBO4SThQAElEv3+CfQ2lcamV+bvPd0N3BiSqkR+g07xiYhIWFIPSkREwpICSkREwpICSkREwpICSkREwpICSkREwtL/A1A3LRm1nBHQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_theta(results):\n",
    "    plot(results.theta, color='C0', label='theta')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angle (rad)')\n",
    "    \n",
    "plot_theta(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot `omega`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd5hU5d3/8fd2emfpXbipu3QLNpAmCNJREGyJxhiTxyS/J1HTTJ5EE6MxJmKNitJUepNiA1SkwwILN733ssD2MvP74wxkRVhmYWbPzOzndV1c7JyZOfM5uux373O+576jvF4vIiIioSba7QAiIiKXogIlIiIhSQVKRERCkgqUiIiEpFi3A/jDGJMAdAEOAwUuxxERkcCJAeoAq6y1OYWfCIsChVOclrkdQkREguYW4KvCG8KlQB0GmDhxIrVr13Y7i4iIBMiRI0cYPXo0+H7OFxYuBaoAoHbt2tSvX9/tLCIiEnjfu3yjJgkREQlJKlAiIhKSVKBERCQklXiBMsZUMcbsM8Y8UNKfLSIi4cONJonXgXqB3KHH4+HAgQNkZGQEcrelXvny5alfvz7R0Rpoi0jJK9ECZYy5H6gEbAzkfk+cOEFUVBTGGP0wDRCPx8PBgwc5ceIEiYmJbscRkVKoxH6aG2OaAL8HHgr0vtPS0qhVq5aKUwBFR0dTq1Ytzpw543YUEQlBHo+HPaf34/F6gvYZJfIT3RgTA0wAfmmtPRLo/RcUFBAXFxfo3ZZ6cXFx5Ofnux1DREJMWtYZ/rTkn/zvor+weEfwJvkpqVN8vwWstXZ6sD4gKioqWLsutfTfVEQulnpsO/9c/h9OZ5+hcplKtKtlgvZZJVWg7gHqGmOG+B5XBMYZY7paa39cQhlEROQqeb1e5tjFTEqZhcfroXXN5vzsxoepWrZy0D6zRAqUtbZl4cfGmPXAy9ba90ri80VE5Oql52bw6orxrDnk9LcNatWHkW0HEBMdE9TPDZe5+ERExAU7Tu7hH8vf5njGScrHleXx6x+gc72kEvlsV9rerLXtS+vo6ZtvvmHo0KF07NiRAQMGsGDBAgB69OjB+PHj6dOnD+3bt+dXv/oVK1eupH///nTo0IFf/OIXFBQ4cymeOXOGp556im7dunHbbbfx0ksvXWhmyMvL449//CNdunShZ8+evPXWWxjz33PEkyZNYsCAAXTq1Ikbb7yRF154oeT/I4hIyPN6vSzY/iW/+/xFjmecpFnVRvy199MlVpwgQkdQzy19lXWHN5XIZ3Wo05anbn3cr9du376dRx99lL/+9a/07t2bVatW8eMf/5iaNWsCMHXqVCZPnkxWVhb9+vVj+/btjB8/npycHAYPHsySJUvo0aMHv/rVr4iPj2fhwoVkZmby05/+lDfeeIPHH3+ccePGsX79eubNm0dMTAw//vF/L/GtXbuWl19+mcmTJ9OsWTNSUlIYNWoUffr0ISmp5L7pRCS0ZeZm8frqCXy7fy0Ava+7lfvbDyMupmS7pXXjUAmaN28e119/Pf369SM2NpYbb7yRAQMGMGPGDADGjh1LtWrVqFevHg0bNmTgwIHUqFGDevXq0bx5cw4cOMCJEyf44osv+N3vfkeFChVITEzk8ccfZ8qUKQDMnj2bRx99lMTERKpXr84TTzxx4fNbtWrFzJkzadasGadPnyY7O5vy5ctz7NgxV/57iEjo2XN6P79e/Bzf7l9L2dgy/M+ND/ODTveWeHGCCB1B+TuiKWmnTp2ibt2639lWv359VqxYAUCVKlUubI+OjqZSpUrfeezxeDh06BAAffv2vfCc1+slLy+PnJwcjh49+p1FHevUqXPh65iYGN544w0WLlxI1apVad26NR5P8G6yE5Hw4fV6WbxzGePXfUyeJ59GVerz85t+SJ2K7s0kE5EFKlTVqVOH1atXf2fb/v37qVGjBrt37/brvqPExESio6NZtmwZZcuWBSA9PZ2TJ0+SkJBAnTp1OHz4MMnJyQAcPXr0wnvfffddUlNTWbRoEZUqVcLr9dKlS5cAHqGIhKPMvCzeXDWRb/avAaBn05t5oMNw4mPjXc2lU3wlqF+/fqxevZr58+dTUFDA8uXLmTNnDgMGDPB7H7Vr16Zr1648//zzZGRkkJ6ezlNPPcVvfvMbAIYMGcKbb77J8ePHOX36NOPGjbvw3nPnzhEXF0dsbCxZWVm89NJLnDt3jtzc3IAfq4iEh12n9vHrRc/xzf41lIlN4Kc3PMQjXUa7XpxABapENWrUiHHjxvH222/TuXNnnn32WZ599lluvvnmYu3nxRdfJD09nV69etGjRw+ioqJ4+eWXAXj44Ydp3bo1ffr0Yfjw4bRp0+bCNFAPPfQQZcuWpVu3bvTs2ZMTJ07QrVs3tm/fHvBjFZHQdr5L7zefvcCR9OM0qlKf53s/xc2NQuesSpTX63U7wxUZYxoDuz/77DPq16//vee3bNlCq1atSjxXKNqwYQONGzemcmXn7u4lS5bwzDPP8NVXX13V/vTfViTyZORm8tqqD1h5YD3gdOmNbT+MeBcaIQ4cOMAdd9wB0MRau6fwc7oGFWGmTZtGVlYWf/7zn8nOzub999/nlltucTuWiISIbSd28c/l/+F45inKxpXhsS5juKFBR7djXZJO8UWYJ598kpycHG655RZ69uxJjRo1ePrpp92OJSIu83g9zNqyiN9//iLHM0/RrFoj/tb76ZAtTqARVMSpWrUqr7zyitsxRCSEpGWf5dUV49lwJBWAu0xPRrW7m9iY0C4BoZ1ORESuScqRLfxrxXucyT5LxfjyPH79/XSs287tWH6JmALl9Xq1flGAhUMDjYhcWr6ngI82zWHWlkV48dImsQVPXP8g1cpVufKbQ0REFKiYmBjy8vKIj3e/bz+S5OXlERsbEd8iIqXK0fTj/HP5O+w4tYeoqChGthnA4FZ9iY4Or7aDiPjpU6VKFY4ePUq9evXC7n9AqPJ4PBw9evRCu7qIhIev9q7irdWTyMrPpnq5qvzshodoWfM6t2NdlYgoUDVq1ODAgQNYa92OElHKly9PjRo13I4hIn7Iysvm3bUf8eWe5QBcX78Dj3YZTYX48i4nu3oRUaCio6Np2LCh2zFERFyx89Re/rn8PxxJP05cTBwPdhjOHU1vDvvr8hFRoERESiOP18Nc+ymTU2ZR4PXQsHI9/ufGh6lfuc6V3xwGVKBERMLQqcw0Xl35HhuPOpc27mzendHJg12ZrihYVKBERMLMqoMbeH3lB5zLzaBSQgUe6zqWTmFyb1NxqECJiISJ7Pwc3l8/jU93LgMguXZrHu86liplI7PbVgVKRCQM7Dq1l1e+fZdD544SGx3L6KRB3NmiO9FRkXtrjQqUiEgI83g8zLGfMmXTbAo8BTSoVIef3vgQjap8f+mhSKMCJSISok5knOLfK94j9bizqGjf5rdzX9LgkFjttiSoQImIhKCv963irdWTyczLonJCRR7rOpaOddu6HatEqUCJiISQzNws/rN2Csv2rgSgY912PNblPiqXqeRyspKnAiUiEiJSj23j3yvGcyLzFPExcYxtP4xezW4J+xkhrpZfBcoYUx24A+gMJAIFwBFgFbDYWpsRtIQiIhEuryCPDzfNZc7WxXjx0qxqI5644QHqVqrtdjRXFVmgjDHXAc8A9wLHgVTgJBAD3AA8AlQwxkwA/mat3R7cuCIikWVf2kH+teI99qYdICoqiqGt+jG0TT9io2Pcjua6yxYoY8zvgLHAeKCTtXbzZV7XBhgNLDbGvGet/UMwgoqIRBKP18P8bV8wOWUmeZ58alWoyRPXP0CLGk3djhYyihpBZQCtrbW5Re3AV7ieNsb8CXg8kOFERCLRiYxTvLpyPJuPbQPgjqY3c3/7oZSJK+NystBy2QJlrX2xODuy1mYBf7/mRCIiEcrr9bJ0zwreWfchWXnZVE6oyKNdRtO5XrLb0UKSv00S0cD9wKfW2v3GmF8D9wErgZ9Za88FMaOISNg7m32ON9dMYuWB9QB0qZfMo51HU6lMRZeThS5/28yfA34A9PY1TvwJ+BvQC/iH7zkREbmE1Qc38MaqiZzJOUfZ2DI82HEEtzW+odS2j/vL3wJ1HzDcWrvGGPMmsNRa+4wxZhqwCBUoEZHvyczLYvy6qXyx+xsA2iS24Mddx1KzfHWXk4UHfwtUVeB8C3k/nFETwBmgdEwKJSJSDJuOWl5b+T7HM08RFx3LqFIw+3ig+VugNgJjjTFHgLrALGNMHPALYH2wwomIhJvc/Fwmpcxk/vYvAGhatSE/uf6BiFmGvST5W6B+CcwAqgF/sdbuMMaMA4YB/YMVTkQknGw7sYtXV47n8LljxERFM6T1nQxufaduur1KRd2oW89aexDAWrvMGJMIVLbWnva95Dng59ba7BLIKSISsvIK8vh48zxmbV2E1+ulfqU6/OT6+2larZHb0cJaUSOolcaYNGAxTiPEl4WKE9ba/cEOJyIS6naf3s+rK8az78xBoohiYMvejGh7F/ExcW5HC3tF3ahbzzeNUW/gJ8BEY8x6fAXLWru6hDKKiIScfE8BM1I/YXrqJxR4PdSuUJPHr78fU6OZ29EiRpHXoHzTGG0G/mGMiQe64RSs140xDYEvcGYzfzvoSUVEQsS+tIO8umI8u9OcE0l9m9/OqKRBlIlNcDlZZPF7PSjfnHxf+P485VuCo6fvjwqUiES8Ak8Bs7Yu4uPN8yjwFFCzfHV+3HUsbRJbuB0tIhXVJHGrH+8/DHzg74cZY+4C/gI0AY7hLNHxhr/vFxFxy4Ezh3l1xXh2nt4LQK9mt3Bf8hDKaoLXoClqBPXlRY+9QBTgwVmwMM73dS5Q7kofZIypA0wFBltrPzHGdAS+NsasstauvYrsIiJBV+ApYPbWxXy8eR75nnyql6vKY13GkFS7ldvRIl5RBarwDIbDgCeBR4HV1toCY0wS8Abwnj8fZK09bIypaa0955t8tjqQD2iiWREJSfvPHGLcyvfZecoZNfVo2o2x7YdSLq6sy8lKh6K6+C4s4+5b62motXZVoedTjDGPA5/gFKor8hWncjhTJMUCf9UqvCISai41avpRl/tIrt3a7Wilir9NEpVwlnm/WMVi7OO8bKA8kATMN8Zst9b+p5j7EBEJin1pBxm38n12nd4H+EZNyUMpF69RU0nzt7hMBd41xjwJrMO5FnU98BLwfnE+0Fp7/rrVat/M6HcDKlAi4qp8TwEztyxkWup8CjwF1ChXjR91uU/Xmlzkb4F6AngNmFXoPXk47eX/688OjDG3AS9ZazsV2pwApPmZQUQkKHaf3s+4le+zN+0AoA69UOFXgfIt5/6AMeYJwPg2b7XWphfjs9YD9YwxPwf+iTMCexgYXIx9iIgETF5BHtNS5zNzyyI8Xg+J5avzoy730bZWS7ejCcW4fuSbLLY5/70W1ckYkwB0stY+d6X3W2vPGGP6Aa8Avwf2Az+w1i4pfmwRkWuz/eRuXlv5AQfOHiaKKO5s3p172w2kjEZNIcOvAmWMeRgYh3Pv0/n7ofB9nYIzs/kV+e53urn4MUVEAiMnP5cPN81h3rbP8Hq91KmQyI+63kerms3djiYX8XcE9RTwb5xZILYAN+GsDfUemuZIRMJE6rFtvLZqAkfTjxMV5Zt5vE1/4mO1MHgo8rdANQBetdaeNMasA9pYa2cZY34GvAy8GrSEIiLXKDMvi4kbZrB45zIAGlSuy2NdxnBd9cbuBpMi+VugzgDnbwLYBiTjdPRZoHHgY4mIBMbaQxt5a/VkTmadJiY6hiGt+jK4VV9iY4p7C6eUNH//Dy0GXjLGPAp8A/zWGDMeGAUcDVY4EZGrdTb7HO+t+5iv9jkT4DSr1ojHuoyhYZV6LicTf/lboJ7EuSH3LuB14IfAbpy59B4JTjQRkeLzer18vW8V7677mHM56cTHxHFPu4H0a96D6Ohot+NJMfhboNrizEKe5XvcwxjTGkiz1h4KTjQRkeI5kXGKt9ZMZt3hTQC0TTQ80mU0tSvUdDmZXA1/C9RHQA+clnIArLWpQUkkIlJMHq+HRTuWMillJtn5OZSLK8vY9kPp3uQmoqKirrwDCUn+FqgdOI0RKVd6oYhISTpw5jBvrJqAPbkLgK712/Nwx3uoWrayy8nkWvlboLYD7xljngJ2AlmFn7TWjgh0MBGRouQV5DFzy0Kmb1lAgaeAqmUq81CnkVxfv4Pb0SRA/C1Q+RRz1nIRkWCxJ3by+qoJHDx7BICeTW9mdPJgysdfcXFvCSP+Thb7YLCDiIhcSWZeFpNSZrJ4xzK8eKlTMZFHO4+mdWILt6NJEFy2QBljPgWettau9GdHxphuwJ+stT0CFU5E5LyVB9bzztoPOZWVRkxUNHe36sOQ1v2Ij4lzO5oESVEjqF8DbxhjCoDpwAIg1VqbC+CbyTwZuA24z/eeHwYxq4iUQqey0nhn7YesPLAegObVGvNIl9E0qlLf5WQSbJctUNba1caYLjjrNf0E+BMQZYzJwJnNvDzOtamvfM9Ns9Z6gx9ZREoDj9fDpzuXMTFlJll52ZSJTWBU0iB6N7tVN9yWEkVeg/Itzz4NmGaMqQh0AGoBHuAIkGKtPRf0lCJSquxLO8ibqyexzdc63rluEg91GkmNctVcTiYlye/ZEn2FaGkQs4hIKZebn8vU1PnM2bqYAq/nQut413rtdcNtKaTpfEUkJKQc2cJbqydxNOMEUUTRu9mtjEoaRLn4sld+s0QkFSgRcdWZ7LO8v34ay/Y6DcMNKtfl0c6jaVGjqcvJxG0qUCLiCo/Xwxe7vmFCygwycjOJi4ljeJv+3GV6Ehsd43Y8CQF+FShjzP3AdDVEiEgg7D9ziLdWT2LriZ0AJNduxQ863UstzTouhfg7gvoNMM4YMxeYAHxirc0PXiwRiUS5+blMS/2E2VsXUeD1UDmhIg90HM5NDTqrCUK+x9+pjpobY24A7gXeBOKNMR8DE621y4IZUEQiw/rDqfxnzWSOZpwAoGezWxiVdDcV4su7nExCVXHazL8FvjXGPImzNtQgYIEx5gTOqOota+2eoKQUkbB1KiuN8eumsnz/GgAaVq7HI51HqQlCrqhYTRLGmBigDzAS6A+cBWYCzYFNxpj/Z619LeApRSTseDweFu1cyuSNs8jKyyYhJp5hbfrT39yhJgjxi79NEt1xTu8NAcrgFKX7gEW+2SYwxvwP8BygAiVSyu06tZc3V09i1+l9AHSs246HO46kZvnqLieTcOLvCGoR8DnwJE43X8YlXrMWmBSoYCISfjJzs5iyaTYLdyzB6/VSvWxVHuw4gi71ktUEIcXmb4EagzMZbF7hjb4ZzftZa2dYa5eiqZBESiWv18s3+1czft1U0rLPEh0VzV3mDoa36U+ZuDJux5Mw5W+Bmgh8Bhy/aHtTnFGT5iIRKaUOnTvKf9ZMYePRrQC0qN6UH3a+V8thyDUrasHCx4BnfQ+jgFRjzMXLaVQA1gUpm4iEsNz8XGZsWcisrYvI9+RTIb48o5MG0b3pTURHaTkMuXZFjaDeAjKAaOAdnDWfzhR63guk44ysRKQUWXd4E++s+fDCPU23N7mR+5IGU6lMRZeTSSQpasHCfOB9AGPMbuBrzR4hUrqdyDjFe+s+ZuVBZ3XbBpXr8sNO99Ky5nUuJ5NIVNQpvr8Bz/o69voD/Y0xl3yttfZ/gxNPREJBvqeAefYzpm6eR05BLgmxCQxv059+LXroniYJmqJO8XUB4gp9fTla5l0kgm0+to2310zm4NkjANzQoCP3tx9G9XJVXU4mka6oU3zdL/X1ecaYWJ3yE4lcp7PO8MH6aXy1bxUAtSvU5OFO95Bcu7XLyaS08HcmiYrAq8BWa+1ffJt3GWMWA09YazODFVBESlaBp4CFO5bw4aY5ZOVlExcTx5BWfRnYshdxMXFX3oFIgPh7H9SrQBvg5ULbxgB/B14EHgtwLhFxwdbjO3h7zRT2nTkIOFMUPdRhBIkVaricTEojfwtUP6CHtTbl/AZr7RJjzKPAAlSgRMJaWvZZJmyYztI9KwCoWb46D3YYQed6SS4nk9LM3wIVhTNJ7KXEByiLiJSw753Oi45lYMveDG7Vh/hY/dMWd/lboOYArxpjHrTWbgIwxrQC/gXMD1Y4EQmeLce3886aD9nrO53XoU4bHuwwgtoVE11OJuLwt0A9ibPERooxJgentTwBZ5bzJ4KUTUSCIC3rDBM2zGDpXt/pvHLVeKDjCDrXTdKM4xJS/F3y/TRwmzGmNdAayAW2WWu3BjOciAROvqeAhdu/5KPNcy+czru7VW/ubtmHBJ3OkxDk94q6xpgywI043XzRQB1jzBFrbVqwwolIYGw+to131n7I/jOHAKc774EOw6ldoabLyUQuz9/7oFrgnM4rjzN7eRQwCvi9MeZWa+0OP/fTC3geZ4n4Y8AL1to3ria4iFzZyczTfLBhOt/sWw1ArfI1eKDjCDrVbedyMpEr83cE9U+cwjT6/E25xphywHjgH8CAK+3AGNMAmAbcD8wCOgELjTF7rLULryK7iFxGXkEe87Z9zrTUT8jJzyEuJo7Bvptt43WzrYQJfwvUrUDXwjNGWGszjTHPAt/4uY/GwCRr7Qzf41XGmC+BboAKlEiArD+8mXfXfcThc8cA6Fq/PWPbDyOxfHWXk4kUj78F6jRQ5RLbqwB5l9j+PdbaZcCy84+NMdWAW4AP/MwgIkU4mn6c8eumsvqQcz993Yq1eLDjCM2dJ2HL3wI1A3jNGHO/tXYdgDGmI84USNOL+6HGmMrAbGAFzuk+EblKOfm5zNiygDlbF5PnyadMbALD2vSnX/PuxMb43QclEnL8/e59Buf60RrffVDgzCAxA/h5cT7Q13AxC0jFuablKc77RcTh9Xr59sBa3l8/jZOZpwG4tdH1jE4eTNWylV1OJ3Lt/L0P6izQyxjTFuc+qCxgi7/de+cZY27FKU6vA09ba7WWlMhV2Jd2kHfXfcTmY9sAaFKlAQ92HEnLms1cTiYSOEWtqFvuEpt3+f585zX+LLdhjGkGzAWesdb+q/hRRSQ9N4OPNs1l0Y6leLweKsaX5552d3NH025ER0e7HU8koIoaQaVz5dVyo3yv8WfN58eBisBzxpjnCm1/1Vr7Kz/eL1JqeTwePt/9NZM3zuZcTjpRUVH0ue42RrYdQIWE8m7HEwmKogrU91bRvRbW2p9TzOtVIgJbj+/k3bUfsjttPwCtazbnwY4jaFSlvsvJRIKrqCXfl1y8zRhTAWcWiFQg3lp7LojZREq1k5mnmbhhxoUl16uXq8qY5KHc2KCjJnWVUsHfqY7icVbTfcS3qQXwV2NMWZxOvDNByidS6uQW5DHXfsqM1AXkFOQSFxPH3S17aVJXKXX8bTP/E3ATzowSi3zbXgDeAV4CHg58NJHSxev1surgBt5fP5VjGScBuL5+B8a0H6pZIKRU8rdAjQDGWGu/McZ4Aay1K40xP0Q32opcs31pBxm//mM2HrUANKhUhwc7jqBtrZYuJxNxj78FKhE4contZ4FLtaOLiB/Sc3xt4zudtvHy8eUY2XYAvZrdQky0P82xIpHL3wK1FPgZ/1091+u7LvVb4KtgBBOJZAWeAhbvXMZHm+aSnptxoW18RNu7qJhQwe14IiHB3wL1M2CBbz2nMsB7ON18BUDv4EQTiUwpR7Ywfv3UC4sHtk00PNBhOA2r1HM5mUho8Xeqo23GmFY4ixS29r1vIjDBn1kkRASOnDvG+xums/rgBgASy1dnbPthdKmXrLZxkUvwt8381zhrOb0b5DwiESczL4vpqQuYv+1z8j35JMQmMKRVX/qbO7R4oEgR/D3Fdw/wf8aY5cAE4CNr7engxRIJfx6Phy/3LGfyxtmcyT4LwG2Nb+DepLupVvZSy6uJSGH+nuJrb4wxOIXqp8ArxphFwCRgprU2K4gZRcLOluPbeW/dx+w+7UxP1KJ6Ux7oMJzrqjd2N5hIGPF7NTNrrQWeBZ41xrTDuTfqDd+fSsGJJxJejmWcZMKG6Xy7fy0A1ctWZXTyILo17KLrTCLFVKzlNn2t5X1xilM/nKXgJwchl0hYyc7LZsaWhcy1n5LnySc+Jo6BLXszsGUvysQmuB1PJCz52yRxF05RGgjkA1OBQdbapUHMJhLyPF4PS/esYFLKTNJ815m6NezM6OTB1ChXzeV0IuHN3xHUh8AcYAzwibU2P3iRRMLD1uM7eG/dx+w6vQ+A66o15oEOw2lRo6nLyUQig78Fqpa1Nj2oSUTCxLH0E0xImXHhOlO1slUYlTSImxt1ITpKq9qKBIq/XXwqTlLqZeZlMXPLQubaz8j3XWe6u2VvBug6k0hQFKtJQqQ0utT9TDc36sropEFUL1fV5XQikUsFSqQIm45uZfz6aexNOwDofiaRkuRvF99rwN+ttTuDnEckJBw6d5QJ66ez+lAKADXKVWN08iBuatBZ9zOJlBB/R1D3An8LZhCRUJCek8HUzfNYuGMJBV4PZWITGNSqD3e1uIN4LbcuUqL8LVBvAi8ZY/4C7AK+M7WRZjSXcJfvKWDRjiV8vHkeGbmZRBFFjyY3cU+7gVQpW9nteCKlkr8F6kGgOs6NupeipT8lLHm9XtYcSuGDDdM5fO4Y4KzPNLb9MBpXre9yOpHSzd8CNSyoKURcsOf0ft5fP41NxywAdSokMqb9EDrVTdJ1JpEQ4O99UEsu95wxpm7g4ogE36msND7cOIcvdy/Hi5fy8eUY3qY/vZvdSmyMGltFQoW/XXwtgBdwVtM9fzovCkgAEv3dj4ibsvNzmGs/ZdbWxeTk5xATHUOf625jWOt+VEgo73Y8EbmIv4XlNaCa7++/AL8BmgI/AH4YnGgigXF+QtcpG2dzKisNgK712jM6eTB1Kia6nE5ELsffAnUD0N1au9IYMxL41lr7d2PMdmAUoKXgJSRtPraN99dPvbBwYJOqDbi//TBaJ7ZwOZmIXIm/BSoKOOL7eivQAfgKmIkzmhIJKYfOHmHChhkXbrTVhK4i4cffArUBGAK8DGwGbgP+BdTDKV4iIeFsTjpTN81j8c6lFHg9JMQmOBO6mp4k6EZbkbDib4H6A33dBY0AABRoSURBVDDLGJMHTAB+Y4z5EmgJzA9ONBH/5RbksWD7F0xPXUBmXhZRUVHc0fRmRra9SzfaioQpf9vMF/o6+WKttYeMMTcBY4FPgFeCGVCkKF6vl+X71zAxZSbHM04CkFy7NWOSh9CwSj2X04nItfC7Pdxau6/Q15uA/w1KIhE/bT2+kw/WT2X7qT0ANKhclzHJQ2lfp7W7wUQkIC5boIwxqwCvPzux1nYNWCKRKzhy7hgTU2ay4sA6ACqXqcTItgPo3uRGYqI165ZIpChqBDW3xFKI+OFcTjrTNs9n4c6lFHgKiI+JY2DLXgw0vSgTV8bteCISYJctUNbaZ0syiMjlOA0QXzIj9RMy8rKIIoruTW5iZNsBVCtXxe14IhIk/k51VORaUNZaXY+SgPN4PSzfv4ZJKbMuNEC0q9WSMclDNdO4SCngb5NEl0u8rwlQFZgS0EQiQOqx7XywYRo7T+0FzjdADCG5dmvNNC5SSvjbZt79UtuNMS8D+QFNJKXaobNHmJAyk9UHNwBQpUwl7mk3kNsb30h0tGaAEClNrnUW8leA1cAvA5BFSrEz2Wf5ePM8Pt35FZ4LM0D04i7TkzKxCW7HExEXXGuB6g0UBCKIlE45+bnM2/YZs7YsIis/+8IMECPa3kVVzQAhUqr52yRxqXuiKgPXAc8X90ONMV2BudZarXVQSnk8Hpbs+ZYPN825sARGxzptGZ08mAaVtQamiPg/gprHdwuUF8gFVllrP/P3w4wxUcDDwN/9TigRZ/3hVCZumM7eMwcBZwmMMclDaFurpcvJRCSU+Nsk8YcAfd6zQH/g/9AyHaXOntP7mbBhBilHtwBQo1w17m13N90addYSGCLyPf6e4nvnMk+dH0kdBKZZa7dcYVevW2t/Z4y53f+IEu5OZJ5iysbZLNuzEi9eysWVZXCrvtzZojvxMXFuxxOREOXvKb5zwE+AlcBy37YuQDecRQsbAE8ZY4ZaaxdcbifW2kPXkFXCTGZuFjO2LGD+ts/J8+QTEx1Dn+tuY2jrO6mYUMHteCIS4vwtUE2B5621zxTeaIz5LdDRWnuXMeYRnFN3ly1QUjrkF+SzaOdSpm2ez7ncDABuatCJe5LupnaFmi6nE5Fw4W+B6gH8/BLbpwBP+b5eAPwjEKEkPDlrM61lcspMjmacAKBVzesYkzyU66o3djeciIQdfwvUfqAPsP2i7X2BI76vGwGnA5RLwkzqsW1M2DCDHb61mepVqs3opEF0qpukqYlE5Kr4W6B+B0zwNTesAqKBTsAA4EFjTGtgIjA5GCEldB04c5iJKTNYc2gj4ExNNEJrM4lIAPjbZv6RMeYA8DhwH5AHbAJustau8d14+zdgnJ/7+xLQOglh7FRWGh9vmsfnu7/G6/VSJjaBgS17c5e5Q1MTiUhAFGfJ92+Aby7z3EqcDj+JcJl5Wczeupi59lNyC/KIjoqm13W3MKxNf6qUqeR2PBGJIP7eB1UW+BHOab044DsXFay1IwIfTUJJfkE+n+76iqmb53E2Jx2ArvXbM6rd3dStVNvldCISifwdQb0JDMHp1DsbvDgSarxeL98eWMvklFkcST8OgKnRjPuSB2NqNHM5nYhEMn8LVF9glLV2VjDDSGi5uDOvbsVajE4eTGd15olICfC3QOUB24IZRELH/jOHmJgyk7WFOvOGt7mLHk1vUmeeiJQYfwvUP4DnjTE/sNYeD2Ygcc/JzNN8tGkuX+5Zrs48EXGdvwVqBJAEHDHGnMOZIPYCresU3jJyM5m1dRHztn1OXkEeMVHR9LzuVoa16UdldeaJiEv8LVD/vsz26lxUrCR85BXksXDHUqanfkK6b868G+p35N6ku6lTUb9ziIi7/L1Rd3zhx8aY3sBDwCDfPi5XwCQEebwevt67mimbZnM84yQArWo2577kwTSv3sTldCIiDr9v1DXGNAYeBB4A6gPpOO3nKk5hZMORVCZumMGetAMANKhUh9HJg+lQp60680QkpBRZoIwxCcAwnNHSbYAH+BKoB9xqrd0Q7IASGLtO7WNiygw2Ht0KQPWyVRnR9i5ua3wD0dFazVZEQs9lC5QxZhxwLxAPLAYeBmZba08bY/JwWs8lxB1NP86UjbP5et9qAMrFlWVQqz70a96d+Nh4l9OJiFxeUSOoH+Hc+/RnYL619mTJRJJAOJt9jmmpn7Bo51IKPAXERsfSt/ntDG7VR6vZikhYKKpAdQdGA/8E3jHGfA1MA2aURDC5Otn5OcyznzF762Ky8rOJIopbG13PyHYDqFm+utvxRET8dtkCZa1dAiwxxvwEZ92n0cALwMu+lwwxxhyw1mpuvhCQ7yng811fM3XzPNKynf8lHeq0YVTSIBpVqe9yOhGR4rtiF5+1Nhdn5DTNGFMFGIlTrJ4Ffm2MmWKt/UFwY8rleL1eVh5cz6SUmRw+dwyAZtUaMTppMG1rGZfTiYhcPb/bzAGstWnAG8AbxpiGwBhgVDCCyZWlHtvGxA0z2O6bzLVOhUTuTbqb6+t3UMu4iIS9YhWowqy1+3AaKP4cuDjij31pB5mUMpO1hzcBULlMJYa17scdzW4mVpO5ikiEuOoCJSXveMZJPtw0h2V7VuKl0GSuLXpQJq6M2/FERAJKBSoMnMtJZ0bqAhbsWEK+J5+Y6Bh6Nb2FoW3u1GSuIhKxVKBCWHZ+DvO3fc6srYvIyssGoFvDzoxsN5DaFWq6nE5EJLhUoEJQgaeAL3Z/w8eb5nE6+wwAybVbcW+7QTSt1tDldCIiJUMFKoR4vV5WHFjH5I2zLrSMN63akFFJg0iq3crldCIiJUsFKkRsPraNSYVaxmtVqMm97QZyQ4OOREdpMlcRKX1UoFy2N+0Ak1Jmsu7wZkAt4yIi56lAueRYxkk+3Dibr/auUsu4iMglqECVsLM56UxP/YRFO5ZeaBnv3exWhra+k0plKrodT0QkZKhAlZBLzTJ+c6Ou3NN2AIkVargdT0Qk5KhABZkzy/hXfLx5Pmd8s4y3r92aUUmDaFy1gcvpRERClwpUkHi8Hr7dv5YpG2dzJP04ANdVa8yopEGaZVxExA8qUEGQcmQLE1NmsPv0fgDqVEzk3naaZVxEpDhUoAJo16m9TEyZycajWwGoWqYyw9v2p3uTm4hRy7iISLGoQAXAkXPHmLJxNt/sXwNAubiy3N2yN/1a9CAhNt7ldCIi4UkF6hqkZZ1h6ub5fLbrKwq8HuKiY+nb/HYGt+pLhYTybscTEQlrKlBXITM3i9l2EfPs5+QU5BIVFUX3JjcxvG1/apSr5nY8EZGIoAJVDHkFeSzcsZQZqZ9wLjcDgM71khnV7m7qV67jcjoRkciiAuUHj8fD0r0r+GjTXE5kngKgVc3rGJU0CFOjmcvpREQikwpUEbxeL2sObWTyxlnsP3MIgAaV6zI6aRAd6rRVy7iISBCpQF3G1uM7mZgyA3tiJwA1y1VjRNsB3NKoK9HRWv5CRCTYVKAusi/tIJM3zmLNoY0AVEyowJBWfel93a3ExcS5nE5EpPRQgfI5nnGSDzfNYdmelXjxkhCbwF0t7mBAy56UiyvrdjwRkVKn1BeosznpzEhdwMIdS5zlL6Ki6dnsVoa26UeVMpXcjiciUmqVaIEyxiQDrwNJwC7gIWvtqpLMcN6F5S/sYrLysgHo1rAzI9sNpHaFmm5EEhGRQkqsQBlj4oFZwMvArcBQYJExppG19mxJ5cj3FPDZzq+Ymvrf5S+SfctfNNHyFyIiIaMkR1C3A3HW2pd9j6cYY34CjATeCvaHe7welu9fw5SNcziq5S9EREJeSRao1sCWi7ZtBdqVxIe/uXoSn+/6GoC6FWtxT7uBWv5CRCSElWSBqgBkXrQtEyhXEh+eEBNPnYqJDDC96N7kRi1/ISIS4kqyQGUAF/drlwPSS+LDH+w4oiQ+RkREAqQkp0RIBS6+2NPSt11EROQ7SnIE9QUQZYx5Evg3ThdfEjCjBDOIiEiYKLERlLU2F7gTpzCdAp4BBllrj5dUBhERCR8leqOutXYTcHNJfqaIiIQnTcstIiIhSQVKRERCkgqUiIiEpHCZzTwG4MiRI27nEBGRACr0c/17syeES4GqAzB69Gi3c4iISHDUAXYW3hAuBWoVcAtwGChwOYuIiARODE5x+t7SS1Fer7fk44iIiFyBmiRERCQkqUCJiEhIUoESEZGQpAIlIiIhSQVKRERCkgqUiIiEJBUoEREJSSpQIiISksJlJolrYoxJBl7HWcF3F/CQtfZ7dy2HE2NML+B5oDlwDHjBWvuGMSYeZ8XiYTizbrxkrX3OvaTXzhhTBUgBfmetfS+SjtEYUwd4DegOZANvWmt/G2HHeAPwCmCA48Dz1tq3I+EYjTFdgbnW2kTf4yKPyRgzAvgLzswJS4AHrLXHSjx4MVziGBOBfwJ3AFHAJ8DPrLWnfc8H7BgjfgTl+4aZBXwIVAH+DCwyxlRyNdg1MMY0AKYB/4dzTPcCzxlj+gDP4vwgaAZ0Ae43xox1K2uAvA7UK/Q4ko5xFs4UXrWAG3COZRQRcozGmGicY3zFWlsZ53v1375fGsP2GI0xUcaYHwCLgPhCT132mIwxrYH/AA8A1YHtwJQSjF0sRRzj20A+0ATnF+SqwKu+9wT0GCO+QAG3A3HW2pettXnW2inAZmCku7GuSWNgkrV2hrXW4xsNfgl0A+4H/mytPW2t3QP8HXjUraDXyhhzP1AJ2Fhoc0QcozHmeqAp8FNrbba1djfO9+sXRMgx4vzwSgSijDFRgBfnh1su4X2MzwKP4fySWFhRx3QfMMda+5W1Nht4CuhmjGleQpmL63vH6PuFwwM8a63NsNamAW/x35XSA3qMpaFAtQa2XLRtK9DOhSwBYa1dZq390fnHxphqOJPprsMZVqcWennYHqsxpgnwe+ChQtuqEDnH2Amn8P7BGHPQGLMTGAxkESHHaK09iXPKazyQhzMh6NM4o8ZwPsbXrbWdgNXnN/jxvdm68HPW2kxgP6F7zN87Rt8vxIOstTsKvW4Qzs8eCPAxloZrUBWAzIu2ZQLlXMgScMaYysBsYAWwxre58PGG5bEaY2KACcAvrbVHjDHnn6rg+zvsjxE4/4vFEpyRVEtgAc51GoiAY/T9xp0NjMI5LX0TMB1I870kLI/RWnvoEpuv9L0ZVj+LLnOM32GM+SVOgbrJtymgx1gaClQGUPaibeWAdBeyBJQxpgXO+f1UYDT/Pc7Cxxuux/pbwFprp1+0PcP3dyQcYw5w1lr7B9/jDcaYt3FOE0FkHOMQoJu19v/5Hi8xxvyHyDrG8670vRkxP4uMMXHAv4ABQA9r7VbfUwE9xtJwii8V56JlYS357jA87BhjbsUZNc0EhvmuYZwGjvDd4w3XY70HGGaMSTPGpOGcIhiH0+QSKce4FSjna+Q5LxaIpP+PDYCEi7bl44wSI+UYAfDj3993fhYZY8oBDQmzYzbGVAQW4zSBdLXWri/0dECPsTSMoL7AuUD7JM658KE47eYzXE11DYwxzYC5wDPW2n9d9PQHwO+NMSk4w+1f4rSEhhVrbcvCj40x64GXfW3m6UTAMeL8Iz8OvGiM+QXOP+yHcS5M7yIyjnERTofpIzgX0zsCPwR+AOwjMo6xsKL+/U0CvjLG3A4sB54D1llrt7kR9BpMwRnc3OK7xlRYQI8x4kdQ1tpc4E6cwnQKeAYYZK09XuQbQ9vjQEWcf/jphf78FfgdsAmnU3EVznn/192LGhQRcYy+LqfbcK4/Hca5/vQ3a+00IucYN+Oc5nsU57rTJODX1tpZRMgxXuSyx2St3YjT8PM6cAJoAwx3J+bVMcYkAf2ArsCxQj97DkDgj1Er6oqISEiK+BGUiIiEJxUoEREJSSpQIiISklSgREQkJKlAiYhISFKBEhGRkFQabtQVCQhjzHv8d4qeS3kWZ1b5L4CK1toSmcLGN2/h18DYom6I9M2L9y0wxlprSyKbyLXQCErEfz/Dma26Ds6yGODcsHh+29+Bb3xfZ1zi/cHyU2DDle7Wt9Z6gD8S/jfDSimhG3VFroIxpi3OUhlNfOv+uJWjDM6UQT2stZv8fM9O4GFr7ZfBzCZyrXSKTySAfHOQXTjFZ4zx4qwi+xTOXHurcRZ1+3/AGOAs8JS19gPf+ysCL+IsGe4FPsdZTvtySx/cA6QVLk7GmN8CjwA1cdZCe9pa+0mh98zAGQ1+GYBDFgkaneITCb7ngf/BWdK9IbAWpzB1wVkb6Q1jzPm1hN7EKWR9cObp8wILjTGX+2WyP84cfgAYYwb7Pus+nJm05wEfG2MqFXrPAqBnEfsUCQkqUCLB96q19gvfsgRzcdbGedrXqPASzvo5TYwxTXFGRKOstat8o6IxQGOg72X23RlnYtLzGuOsM7XXd+rxjziTteYVek0qzkzb35kxXiTU6DcokeArvDx2JrDHWnv+4m+27+8EoJHva1toBWFwFnwzOMXtYrVwZo0+bwJOp+EuY8wanNWW37XWZhV6zUnf34nFPA6REqURlEjw5V302HOZ18X6XtsBaF/oTwvg3cu8xwNEnX/gW0amE86I6xvgASDF19Rx3vl/9wV+H4GIC1SgRELHFiAOKG+t3WGt3YGzTtQLOEXqUo7gNEMAYIwZAjxqrV1krf0ZzsjrHM4aPufVLPRekZClU3wiIcJaa40xs4H3jTGP46y2+2ec5oqtl3nbGiC50OMY4AVjzFGcjsEbgNq+r89LxllWvvCpR5GQoxGUSGi5H6eYzMRZkbUy0Mtam3aZ18/D6fYDwFr7MfB7nFHXNuD/gJ9Yaz8v9J5bgQXWWp3ik5CmG3VFwpgxphywB+hrrV3rx+ujgb04nYLLghxP5JpoBCUSxqy1mTijpcf9fMvdwC4VJwkHKlAi4e8fQJK5qDf9Yr7R0zPAj0oklcg10ik+EREJSRpBiYhISFKBEhGRkKQCJSIiIUkFSkREQpIKlIiIhKT/D6UNrCVNdahuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_omega(results):\n",
    "    plot(results.omega, color='C2', label='omega')\n",
    "\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angular velocity (rad/s)')\n",
    "    \n",
    "plot_omega(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot `y`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhU1f3H8fedyWSyJ0D2BGTTg6K44L5X3FdQ3K3i3tb2Z0tdqpaiVlyrtbXu2mLdtVXrilr3XUTBlaMoKCEJSQjZ92R+f9wBIwWcQDJ3ZvJ5Pc88ydw7d+ZzHnG+Ofeee44TCoUQERGJNT6vA4iIiKyNCpSIiMQkFSgREYlJKlAiIhKTkrwOsLGMMUFgB6AC6PY4joiI9I0fKALmWmvbe++I+wKFW5ze8DqEiIhslD2AN3tvSIQCVQFw//33U1hY6HUWERHpg8rKSk488UQIf5f3lggFqhugsLCQ0tJSr7OIiMiG+Z9LNBokISIiMUkFSkREYpIKlIiIxCQVKBERiUmJMEhCRCRh9PT0UFZWRnNzs9dR+kUgECA/P5+srKw+H6sCJSISQ2pqanAcB2MMPl98n+QKhUK0traybNkygD4XKRUooKe9hbYy6z5xHHAcHJ8fx+cHnx/HH8BJWvVIxhdIwQkku/tFRPpRXV0dI0eOjPviBOA4DmlpaZSUlFBeXq4CtSGq/vMXWr76oM/HOUnJ+IKpOMmp+IJp+FLS8aek40vJwJeagT8ty32kZuFPz8afMQR/ehaOPzAArRCRRNDd3U0gkFjfEampqXR2dvb5OBUoIHPCPoR6uiHUAyEIhXqgp9vd1t1NqLvTfXR10tPZTmjVo6uD7q4OaK7v0+f50rJIyhiKP3MoSZlDScrOIykrl6TsXJKy80nKGqbemcgg5jiO1xH61Ya2RwUKSB+3E+njdurTMaFQD6HODno6Wulpb6WnvYWetmZ62proaW2iu7WRntZGulsa3EdzPd3NdXQ319PT0kBHSwNULVn7m/v8JGUNIzCkkKQhhQSGFBEYWkRgWDGBnAIcv/6ziUji0zfdBnIcH05yCr7kFMgYEvFxoZ5ut2A11tLVWEtXQ833j/pquuqq6W6qpauuiq66Klj88Q/fwOcnkFNAILeU5NxSkvNGEMgbTnJuiU4dikhCUYGKMsfnJyljCEkZQwgWjVnra3q6OtxitbKSzpWVdNZW0llbTueKcrrqq93fa8tp+fL97w/y+UnOLSE5fyTJBSMJFowiuWAU/rTMKLVMRKR/qUDFIF9SMsnDSkgeVvI/+3o62+msraCzZikd1WV01Cyls/o7Omsr6aj6jo6q7+DT11e/Pik7j2DRGIJFY92fxWPxBdOi2RwRSUAzZsygra2N6667bvW2/fbbj4suuoh99tmnXz5DBSrO+AJBggUjCRaM/MH2no42OqqX0lG1hI7KxbQvX0zH8iVuT6y+muaF74Zf6RDILSFYvBkppYaU4eMIDCvGceJ/SKtIIqp4aBatX38Ytc9LHbMdRcdd8qOvO+KIIzjrrLNob28nGAyyYMECGhsb2WOPPfotiwpUgvAlp5BSsikpJZuu3hbq6aazZhntlV/TXr6I9oqvaa9cTGdNGZ01ZTR9/LJ7bEqGW6xGbEHKiC0IFo7WQAwRWa+JEyeSk5PDq6++ygEHHMDTTz/NwQcf3K9D5PUtlMAcn5/k/BEk548gc8JPAPf6VsfyJbQv+5K2soW0LbV0N9XSsmgeLYvmuccFgqQMH0fqJluSssmWBIvGaNi7iEci6c14wXEcDjvsMJ555hn2228/nnvuOW6++eZ+/QwVqEHGl5RMSslmpJRsRvaOhxIKheiqr6Zt6Re0ffc5bUs/p3NFOa3fLKD1mwUAOME0UjcZT+rICaSOmkBgWEnC3achIn13xBFHcOSRR/L222+TlpbG1ltv3a/vrwI1yDmOQyAnn0BOPplb7QVAV9NK2r79jNZvP6V1ySd0rayk5cu5tHw5F4CkrFxSR29D2phtSR01QYMuRAap0aNHM2bMGK677joOP/zwfn9/FSj5H0kZQ8gYvzsZ43cHoLO+itbFn9C65GNaF39MV0MNjfP/S+P8/4LPT8rwzUkbux1pYyeqdyUyyBxxxBHMmjWLm266qd/fWwVKflQgO5/ANpPI2mYSoVAPHRXf0PLNfFq/mU9bmaXt209p+/ZTal/6J4GhRaRtuj1pm+1ASuk4XbsSSXBFRUVst912jBgxot/fWwVK+sRxfASLxxIsHsuQ3afS3dpI6+KPafnqA1q+/pDO2grq33uK+veewpeWRdrY7Uk3O5I6emt8SclexxeRftLY2MiyZcu44447OP744wfkM1SgZKP4UzPJ2GI3MrbYjVBPN21llpav5tJs36drZSVNH79M08cv4ySnkrbpRDLG7ULqmG3xBYJeRxeRjbB48WJOPvlk9txzT4444ogB+YyoFyhjTA7wMfAHa+3s8LaLgXOADOBd4Axr7dJoZ5ON4/j8pI7YgtQRWzB0n5PprFlKs32f5oXv0rF8Mc2fvUnzZ2/iBFJI22x7MjbfjdQx26hnJRKHJkyYwPz58wf0M7zoQd0GrJ7DxxhzDnAqsCdQBtwM3A3s70E26SeO45CcN4LkvBEM2X0qnSsraV74Ls1fvEN7xaLVxcoXTCN93M5kjN+DlE3G65qViKwW1QJljDkFyAI+6bX5l8AF1tqvw6+ZDmwSzVwy8AJDCsnZZTI5u0x2i9UX79D0+Vt0LF9M44KXaVzwMv70HDK23IOMrfb+n6mcRAaTUCiUUKNhQ6HQBh0XtQJljBkFzAR2BeaEt6UD44AcY8xHuD2rV3GLliSowJBCcnadQs6uU+ioKaPpszdp+uwNulZWrh5gkZw/goytfkLGlnuQ1IflTETind/vp7Ozk+TkxDn13draukFTIEVlhlBjjB+4DzjPWlvZa9eqb54zgcOBsUAo/FoZBJJzSxm613EM//nfKJ52FVkTD8SXmkFH1XfUvnQP3/31LCofvpLmhe8S6u77ktEi8SYnJ4fly5fT09PjdZSNFgqFaGlpYdmyZeTn5/f5+Gj1oGYA1lr72Brb28M/r141KCI8YOIrY0ymtbYxSvnEY47jrJ6Cadh+02hZ9CGNH79Cy6IPV88T6EvLInOrvcncZhLJuaVeRxYZELm5uZSVlWGt9TpKvwgEAhQUFJCVldXnY6NVoI4Dio0xR4afZwK3ALOBFXzfk+qdKXFOwEqfOP4A6WYn0s1OdDfX0/jp6zQueInO6qXUv/ck9e89Scrwzcncdj/SN99FowAlofh8vgG56TUeRaVAWWvH9X5ujJkP3GitnW2MaQYuMca8CtQAs4BnrLUN0cgmsc2fnk3OToeRveOhtJd/ReP8l2j6/E13ctulX7Dixb+TOeEnZG23P4GhxV7HFZF+FAs36l4MtAGvAUOBF4HTPE0kMecHpwD3nUbTZ2/Q8NGLdFR+s3pgReqorcmaeABpm26v4eoiCcDZ0OF/scIYMxJY/NJLL1FaqusSg017+SIaPnyeps/eJNTVAbjL3GdNPJDMbSbhT830OKGIrE9ZWRmTJk0CGGWtXdJ7Xyz0oEQ2WLB4LHnFYxk66RQaP36Fhnlz6FpZSe3L97Ly9YfJ2Govsnc4hOS84V5HFZE+UoGShOBPzQhfqzqE1kUfUf/Bs7R+M5/Gj16k8aMXSR29Ndk7Hkbq6G0S6gZIkUSmAiUJxXF8pG06kbRNJ9JRU0bD3Gdp/OTV1SsEB/KGk73jYWRuuSdOUt9vHBSR6InKjboiXkjOLSX3oLMY8avbGfqTk/BnDKWzeik1z9zCd3/7GXVvP0Z3W7PXMUVkHdSDkoTnT80kZ9cpZO90KE2fv0X9u0/RUbWE2lfuZ+Vbj5G17X5k73goSVnDvI4qIr2oQMmg4fgDZG61Nxlb7kXr4gXUvfMEbUs+cW/+/eBZMrfam5xdJhMYWuR1VBFBBUoGIcdxSBu9DWmjt6G9fBF17zxB88J3aZz/XxoXvEzGFruRs9uRJOfpbn4RL6lAyaAWLB5LwVHn0bFiGXVvP0HTp6/R9NkbNH32BunjdiZnt6kEC0d5HVNkUFKBEgGSh5WQf9g5DN3zGOreeYLG+S+5CywufJe0zXZkyB5HEywc7XVMkUFFBUqkl6TsPHIPPJOc3Y6i/t3/0PDhC7R8+T4tX75P2mY7MGSPY9WjEokSFSiRtUjKHMqw/U4le5fJbqGa9zwtX86l5cu5pI/bmSF7HqtrVCIDTAVKZD2SMoYwbN9pZO98hHvqb97z4VN/75ExfneG7HmsRv2JDBDdqCsSgaSMIeTudyrDf3EzWRMPBJ+fps/eYOnt51L97O10NazwOqJIwlGBEumDpKxh5B54JsN/fhMZE/aBUIjGj15g6a2/ZMXL99Ld2uR1RJGEoQIlsgECOfnkH3YOpWf9mfRxuxDq6qD+nSdYessvqHv7cXo6272OKBL3VKBENkJybikFR51H8anXkDJyK3ramql95T6W3vorGj9+hVBPt9cRReKWCpRIP0gpHkvRCTMpPH4GyQWj6G5cQfVTf2PZ3RfQ8s0Cr+OJxCWN4hPpJ6umUEodNYGmT1+n9tUH6ahaQuWDl5M6ZluGTTpFCyeK9IEKlEg/cxwfmVvtTfq4XWiY+ywr3/o3rV9/RNk3C8jadj+G7Hks/vRsr2OKxDyd4hMZIL5AkJxdpzDiFzeTtd0BADR8+DxLb/0lde89Rai70+OEIrFNBUpkgPnTs8k96CxKz7yB1NHb0NPeQu1/Z1N2x3RaFs3zOp5IzFKBEomS5LzhFB73ewqPvZjAsGI6a8upfPhKKh++ks7acq/jicQcXYMSiSLHcUgbO5HUUROo/2AOK994hJZF82j5ZgHZOx3KkN2n4ktO9TqmSExQD0rEA44/QM5OhzH8ZzeRufU+0NPl3uh727k0ff4WoVDI64ginlOBEvFQUkYOeYeeQ/G0qwkWjaG7cQVVj99AxQOX0VFT5nU8EU9FvUAZY3KMMd8ZY6atZd9NxphXo51JxGspJZtSPO0qcg86G19qBm1LPqHszt9S+8r9mjZJBi0velC3ASVrbjTGHAT8PPpxRGKD4/OTtd3+DP/Z38jcZl/o6aLu7ccou/1cmr/6wOt4IlEX1QJljDkFyAI+WWN7HnATcGs084jEIn9aJnmH/JziU64kOX8kXfXVLH/kKir/da2W9ZBBJWoFyhgzCpgJnLaW3X8Hrga+jlYekViXUmooOf1ahu47DSeQQot9j6W3n0v93Gc1Ca0MClEpUMYYP3AfcJ61tnKNfb8Aeqy1d0Uji0g8cXz+8Gi/v5C22Y6EOlpZ8cLdlN9zCe3Ll3gdT2RARasHNQOw1trHem80xmwOnA+cEaUcInEpKSuXwqMvpGDqBfgzh9Je/hXL/n6BBlFIQovWjbrHAcXGmCPDzzOBW4BUoAX4yhgDEAQCxpg6a21OlLKJxI10sxOpI7ei9pX7aZj3PHVvP0bzwnfJPeTnpI7Ywut4Iv0qKgXKWjuu93NjzHzgRmvt7DW2/xqYbK3dOxq5ROKRL5hG7oFnkrHlHlQ/cyudNWVU3DuDrO0OYOg+J+ELpnkdUaRf6EZdkTiVUjqO0tP/RM7uU8Hnd2dKv+M3tHz9kdfRRPqFE+9TqhhjRgKLX3rpJUpLS72OI+KJjqpvqX76Ztor3IGwmVvvw9B9p+FPSfc4mcj6lZWVMWnSJIBR1tolvfepByWSAJLzN6F42lUM/clJOP4AjQtepuz2X9Oy6EOvo4lsMBUokQTh+Pzk7DqFkjP+RLBkM7qbaql8eBbVT99CT1uz1/FE+kwFSiTBJOeWUnzyFQzd56fh3tRLLL1zOi2LF3gdTaRPVKBEEpDj85Ozy2RKTr+OYNFYuhtqqHzgcmrm3ElPR5vX8UQiogIlksCS84ZTPO1Khux1vDvSb94cyu76LW3LvvQ6msiPUoESSXCOz8+Q3adScurVBPJG0LWykvJ7LqH2tQcJdXd5HU9knVSgRAaJYOFoSk+7luydD4dQiLo3/8Wy2RdrYUSJWSpQIoOIkxRg2KRTKDrpMpKycumo/Jpld59Pw7zntcy8xBwVKJFBKHWT8ZSeeQMZW+1NqKuDmjl3sPzRq+lurvc6mshqKlAig5QvJZ38w39F/pTp+FLSafnqA8runK6pkiRmqECJDHIZW+xG6Zk3kLLJeLqb66h86ApW/Hc2oa5Or6PJINen2cyNMWlAPtANVFpr9S9YJAEkZeVSdMJM6t55gpWvPUT9e0/RuuRT8qf8huRhJV7Hk0HqRwtUeKn2XwIHAWaNfZ8ATwN3Wmu/HZCEIhIVjs/PkN2OInXkVlQ98Wc6li9m2d0XkHvgmWRO2NvreDIIrfMUnzFmqDHm78B8YARwPbA7sDmwJbA3cBewGfCJMeYfxpjcAU8sIgMqpWQzSs+4nozxexDqbKP6qZuo+s9f6Glv9TqaDDLr60G9CtwG/MJau665Ud4A/maMyQJODx+zZX8GFJHo8wXTyDviXFJHTaDm+bto+vR12pZ9ScGRvyVYONrreDJIrG+QxM7W2lvWU5xWs9Y2WGv/DOzYf9FExEuO45C59T6UnHYtyfmb0LWykmWzL6L+gzm6Z0qiYp0Fylrb0tc325BjRCS2JeeWUjztKjK32x+6u1jx/J1UPXa9lvCQARfRKD5jzDbAjcB4ILjmfmttVj/nEpEY4gsEyTvobFI32ZLqZ26leeE7tC9frFN+MqAiHWZ+D1AHnAdorn6RQSpji90IFo5m+WPX07F8MeWzL2bY/qeRue1+OI7jdTxJMJEWqLHA9tbaLwYyjIjEvsDQIopPmcWKF/9B40cvUvPc7bQt/YLcg87Gl5zidTxJIJHOJPEGMGEgg4hI/PAFguQd/DPyjjgXJxCk6dPXWTb7d3SsWOZ1NEkgkfagzgLeNcYcAnwD9PTeaa29vL+DiUjsy9xyT4IFI1n+7+vorF7Ksr9fQN6h55Cx+a5eR5MEEGkP6lLcKY4mAocAh/V6HDogyUQkLiTnjaDk1GtJ33xXQh1tVD12vTuXX0+319EkzkXagzoGmGytfXogw4hIfPIFU8mfMp2G4eNY8d97qH/vKdorviF/ynSSMnK8jidxKtIeVC3uqT0RkbVyHIfsHQ6h+KTL8Kfn0PbdZyy7+3zaln3pdTSJU5H2oM4HbjLGnA98DfxgFnPdoCsiq6QM35yS0/9E1ePX07b0C8r/OYPcA04na7v9vY4mcSbSAvU3IAeYu479/kg/0BiTA3wM/MFaO9sYkw/8BZgEOMBzwLnW2pWRvqeIxJakzCEUnXgpK166h4a5z1Lz3O20V3xN7gFn4CQFvI4ncSLSAjW1Hz/zNqD3AjN3AfXAKCAA3AvcDJzQj58pIlHm+JPI3f90gkVjqHn2dhrn/5eO6u8oOOp8kjKHeh1P4sA6C5QxJsdaWwdgrX0tkjczxgxZX8/HGHMKkAV8En7uwx2yfpm1tjm87U7cHpuIJIDMrfYmOXcEy/91De3LvmTZ3y+gYOoFpJRs5nU0iXHrGyTxmjHmQmNM9o+9iTEm1xjze+D19bxmFDATOG3VNmttj7V2srV2Ua+XTgY++vHoIhIvgkWjKTntWlJGbEF300rK751Bw/yXvI4lMW59p/h2A64AyowxbwFzgM+AGtxrRXnA1sBewB648/XttrY3Msb4gfuA86y1lcaYtb0MY8x5uAVKd/mJJBh/ejZFJ8xkxYv/oGHeHGqeuYWOqiUM23caji/iy9gyiKyzQFlrm4BfG2OuBs7GvSa0Ld8PiOjE7ek8A5xprS1fz+fMcN/SPra2ncaYAHAT7o2/+1hrF/a1ISIS+xx/ErkHnkly4ShqnruThrnP0llTRv6U6fhTM72OJzHG6cvCY+FrRsOAHmvtij4ctxAo5vspkjKBdmA2cCHwVHjb4dbaPk3mZYwZCSx+6aWXKC0t7cuhIuKhtqULWf7va+luridpSCGFR/+O5LzhXseSKCsrK2PSpEkAo6y1S3rvi3QUH+BeMwKq+xrAWjuu93NjzHzgxvAw82dwr4XtofupRAaPlOHjKDn1GiofvYaO5YtZds/FFEz+DWljt/M6msSISGeSGBDGmAnAwbhLxVcZY5rCjzIvc4lIdCRl51F8yizSx+1CqL2Fykeuov79p7WkvAB97EH1F2vtNr2eapUzkUHMFwiSf+R0Vr7+CHVvPsqKF/9BR/VScg88E8fvyVeUxAhPe1AiIgCO42PoXseRP2U6TlIyjfP/S8VDV9Dd2uR1NPGQCpSIxIyMLXaj6KTL3clml3xC+eyL6Kxd3wBhSWQR9Z+NMeOBO3Dve0pdc7+1VjcxiEi/SCnZlJJTr6bykavoqPqWZbMvomDqBaSOGO91NImySHtQd4dfezZw+FoeIiL9Jik7j+KTZ5E2diI9rU1U3H85jZ+86nUsibJIr0BuDWxnrf1iIMOIiKziC6ZScPSFrPjvbBrmPkv1kzfRWVvJkD2PxXE0tmowiLQH9TWg6YdFJKocn5/c/U9n2P6ng+Oj7s1HqX7yr4S6On/8YIl765vNfIteT/8O3G2MuRB3Zd3u3q+11n4+MPFERCB7h4MJ5BSw/PEbaPr0dboaV1Bw1AX4UzO8jiYDaH2n+D4FQvzwPqXHe/2+al+IPixYKCKyIdI2nUjxyX+k8uErafv2M8rvuZjC4y4hkFPgdTQZIOs7xTcKGB3+ubbH6F4/RUQGXLBwNCXTriKQN4LOFcson30x7eWLfvxAiUvrLFDW2m9XPXDXcartvS28vQG4IVphRUSSsvMoOfkKUkduRXdzHeX3zaRl0TyvY8kAWN81qK35fmn2U4AXjTH1a7xsS+CAAcomIrJWvpR0Co+7hOpnbqXpk9eofORqcg86i6xt9/M6mvSj9V2Dygae7vX8/rW8pgm4tl8TiYhEwPEHyDvsVyRl5VL31r+pefY2uhpWaBh6AlnfgoWvEz4FaIxZDOxgra2JVjARkR/jOA5D9z6BpKxcaubcSd2bj9LdWEvuwWdrld4EENGNutbaUQMdRERkQ2Vttz/+jCFUPX4DjQteoru5jvwp0/Elp3gdTTZCpHPxLcYdTr6mENABLAMetNbe3Y/ZREQilr7ZDhSdeCmVj1xJy6J5VDxwGYXHXIw/TUvJx6tIZ5L4K5AHPAz8Ovx4AMgFngXmAJcbY6YPREgRkUiklBqKT55FUlYu7cu+pPyfl9BV3+dFwCVGRFqgfgr8zFp7kbX2yfDj98CZwN7W2j8BpwLnDFRQEZFIJOeWUnzKlavvlVp2z8V0VH3ndSzZAJEWqHHAB2vZvgBYNSXSQqCoP0KJiGyMpKxhFJ98BSnDN6e7sZbye2fQVrbQ61jSR5EWqHnA+caY1deswr+fj1ukAHYH9GeKiMQEf0o6hcfPIG2zHehpa6Li/stoWfSh17GkDyItUL8EDgW+M8Y8Z4x5HrcYHQz83BizNzAbuGogQoqIbAhfIEjBUeeTMWEfQl0dVD56NY2fvu51LIlQRAXKWrsA2Ay4FFgEfAb8Hhhjrf0IWAJsb629Z2BiiohsGMfnJ+/QX5C98xHQ0031f/5C/dxnvY4lEYh0wUKstQ24y76vbd+S/gokItLfHMdh2KST8adlUfvyvax44W562lvI2e0ozToRwyK9D2pT4E/ARCDAD5fgwFqb3//RRET6V84uk/GlZFDz3O2sfO1BelobGbrvKThOpFc7JJoi7UHdjjtC7xrcGcxFROJS1rb74ktJp+qJG6l//2m621rIO+RnmhopBkVaoHbAvd9Jc9qLSNzL2HwXfMFUlv/rWpo+fplQRyv5k8/F8Qe8jia9RNqvXQr0y6RWxpgcY8x3xphp4efJxpg7jDG1xphqY8xF/fE5IiLrkzZ6G4qO/wNOMI3mhe9Q+cg19HS2ex1Leom0BzUDuNUY80fgK9z591az1n7eh8+8je/XmQK4DDDAGNwlPuYYY5ZZa//Zh/cUEemzlOHjKD7pMioe/COt33xE5YN/pPDYi/EF07yOJkTeg3oUd3HCh4EPgU+BT3r9jIgx5hQga41jTgFmWWtXhkcD/gk4O9L3FBHZGMHC0RT/9I/4M4fStvQLKh64nO7WRq9jCZEXqFFreYzu9fNHGWNG4S4df1qvbTm4gy9698AWAltFmEtEZKMl55ZSfPIVJOXk017+FRX3zaSrqc7rWINepDfqfmut/RZYAQwFKoHaXtvXyxjjB+4DzrPWVvbalRH+2dJrWwug/rWIRFUgp4Din15BYGgxHVXfUnHfDLoaVngda1CLqECFBzLcAtQBc3GvIf3dGPO0MSY7greYAVhr7WNrbG8O/0zttS0Ndyl5EZGoSsoaRtFP/0hy/iZ0riin/N4ZdNZXeR1r0Ir0FN8fgV2BPYG28LbrgJHADREcfxww1RhTZ4ypwz2FdwswC7c3Znq9dhw/POUnIhI1SRk5FJ10GcmFY+iqW07FP2fQubLyxw+UfhdpgToG+KW19m3CK+taa9/HXQ/qsB872Fo7zlqbZa3Nsdbm4A6S+IW19hfAvcBMY0yuMWYkcF54m4iIJ/ypmRSfOJNgiaGroYbyf86go6bM61iDTqQFKh+3p7OmBjb+etEfcEcDfoZ7+vDfuEPRRUQ840tJp+j4GaSMGE93Uy0V982ko1orCkVTpPdBvQ6cC/wq/DxkjEnGvbb0Zl8/1Fq7Ta/f23BX4tVqvCISU3zBVAqPu4Tlj15D6+IFlN83k6ITZhIsGOl1tEEh0h7U/wGHGGMW4s4oMRt3iY3dgd8MSDIRkRjgCwQpOOZ3pI7elp6WBirun0l75TdexxoUIh1m/hWwOXA1cCPuzbqXAptZa78YsHQiIjHAl5RM4dEXkjZ2Ij2tTVTcfyntFV97HSvh9WU9qHbcntNqxphiY8xUTUskIonOSQpQMPV8lj92PS1fzqXigcsoOv4PBIvHeh0tYW3sIihbAf/ojyAiIrHO8QcoOPK3pJmd6GlrpuKBy2grX+R1rISlVbpERPrA8QcomDKd9HE709Pe4hapZV95HSshqUCJiPSR408if/JvSB+3C6H2FiofvFw9qQGgAiUisml//I0AABKPSURBVAHcIvXr1T2pygcuo11Fql+tc5CEMebgCI6f2I9ZRETiyqqeVNUTf6Z54btUPHi5e59U0RivoyWE9Y3iezrC9wj1RxARkXi0qkgtf/wGWux7VDxwOUUnziRYGNFKRLIe6yxQ1lqd/hMRiYDjT6JgyvTwEPT33SHoJ16mGSc2koqQiEg/cPxJFBw5/fubeR+4jI4qzd23MVSgRET6ieMPUHDU+aSOCU+L9MClmgV9I6hAiYj0I3fGiQtIHbU13c31VNx/KZ21FV7HiksqUCIi/cyXlEzB0ReSssl4uptWUn7/pXTWaWXevlKBEhEZAL5AkMJjLiJYauhuqKHivpl0NazwOlZcUYESERkgvuRUio69hGDRWLrqq6i4/1K6mlZ6HStuqECJiAwgX0o6hcfPILlgFJ215VQ8cDndLY1ex4oLKlAiIgPMn5pB0fEzCOSW0ln9HRUPXk53W7PXsWKeCpSISBT407MpOuFSkoYU0lH5DZUPzaKno9XrWDFNBUpEJEqSModQfOKlJGXl0r7MsvzRa+jp6vA6VsxSgRIRiaKk7DyKTpyJPz2H1iWfUPXY9YS6u7yOFZNUoEREoiwwtJiiE2biS82g5asPqHryr4R6ur2OFXNUoEREPJCcP4LC42bgJKfS/Plb1My5k1BIi0P0pgIlIuKRlOKxFB5zEU5SMo0fvUjty/eqSPWiAiUi4qHUTcZTcOR54PNT/+5/qHv7Ma8jxYz1LVjY74wxhwJXAqOAKuBaa+3txpgM4GZg1Sq+c4BzrLUN0cwnIuKFtE0nkn/4/1H1xI2sfPUBfMF0src/0OtYnotaD8oYUwT8C7jQWpsJHA3caIzZDrgUyMEtXGOB0vA2EZFBIWP87uQedBYAK56/i6bP3vA4kfeiVqCstRVAnrX2OWOMDxgGdAGNgAm/zAn/7AF0B5uIDCpZ2+3P0J+cCISoevImWhZ96HUkT0X1GpS1ttEYkwa0Ay8AN1trvwJuBPYE6oCVQApwRTSziYjEguxdppC98+HQ083yf19H29KFXkfyjBeDJNqAdGAH4DRjzOlAALgHyAMKgSbgDg+yiYh4ynEchu5zMplbTyLU1UHlI1fSUfWt17E8EfUCZa3tsdZ2WGs/wC1CU4EHgVustbXW2irgN8CJxpisaOcTEfGa4zjkHnw2aWYnetqaqXjwj3TWLfc6VtRFc5DEXsaYeWtsDuJeh8oJ/75KFxAK/xQRGXQcn5/8yb8mZZMt6W5aScUDl9PVVOd1rKiKZg9qPlBijJlujPEbY3YFTgf+BLwNXGuMyTbGZAPXAE9aa1uimE9EJKb4kpIpPPpCkgtH07WyksqHrqCnffB8LUZzFF897n1ORwK1uKf3zrDWvoY75LwB+AqwuAMlTotWNhGRWOULplF03O8JDC2iY/liKh+9hlBXp9exoiKqN+paaz8Edl/L9nLg2GhmERGJF/70bAqPn0H57Itp+/ZTqp78C/mTf4Pj83sdbUBpqiMRkTgQyCmg8PgZOME0mr94hxUv/D3h5+1TgRIRiRPBgpEUHv07HH+AhnlzqHvr315HGlAqUCIicSR1k/HkTT4XcFj52oM0LnjZ60gDRgVKRCTOZIzbhWEHnAFA9TO30vLVmnfwJAYVKBGROJS9/YHk7HYUhHpY/vj1tC370utI/U4FSkQkTg3Z63gyJuxDqLOdykeuorO23OtI/UoFSkQkTjmOQ97BZ5M6Zlt6WhqoeGgW3c31XsfqNypQIiJxzPEnUXDkb0kuHOPONvHwlfR0tHkdq1+oQImIxDlfciqFx15MUk4+7RWLqHr8BkI93V7H2mgqUCIiCSApI4fC42bgS82kZdE8aubcGfc38qpAiYgkiORhxRQe8zucpGQaP3qRurcf9zrSRlGBEhFJICml48g74v8Ah5Wv3k/jp697HWmDqUCJiCSYjHG7MGy/aQBUP3UzrUs+8TbQBlKBEhFJQNk7HkrWjodCTxfL/3UtHdVLvY7UZypQIiIJatikk91l49tbqHx4Fl1NK72O1CcqUCIiCcrx+ck/4lyCxZvSVV/N8keuiqt7pFSgREQSmC8QpPCYi8L3SH1N1X9ujJt7pFSgREQSnD89m8Ljfo8vJYOWL+dS+9I/vY4UERUoEZFBIHlYCQVTzwdfEvXvP039B3O8jvSjVKBERAaJ1E22JO+QnwOw4oW7aVkU2+tIqUCJiAwimRP2Jmf3qeF1pG6gffkSryOtkwqUiMggM2TP48gYvwehjjYqH7mKrsbYHH6uAiUiMsg4jkPuob8gWGrobqhh+aNX09PZ7nWs/6ECJSIyCPmSkimceuHqJTqqn7yJUKjH61g/oAIlIjJI+dOzKTzmYpxgGs0L32Hlqw96HekHVKBERAax5LzhFBx5Hjg+6t5+jMaPX/U60mpJ0fwwY8yhwJXAKKAKuNZae3t438XAOUAG8C5whrU2/mY3FBGJM2mjtyb3gNOpmXMn1c/eSmBIASnDN/c6VvR6UMaYIuBfwIXW2kzgaOBGY8x2xphzgFOBPYF8YClwd7SyiYgMdlkTDyRr+4Ohu4vKf11L58pKryNFr0BZayuAPGvtc8YYHzAM6AIagV8CF1hrv7bWtgPTgd9GK5uIiMCw/aaROnobeloaqHzkKnraWzzNE9VrUNbaRmNMGtAOvADcDJQD44AcY8xHxpgq4C5geTSziYgMdo7PT8GU6QRyS+msKWP543/2dGJZLwZJtAHpwA7AacB54e1nAocDY4EQcJ8H2UREBjVfSjqFx1yELzWD1q8/pPYV776Ko16grLU91toOa+0HwB3AzuFdV1trl1prG4CLgX2NMZnRziciMtgFhhRScNT54PNT/+6TNC542ZMc0RwksZcxZs2ZCYNANbACGNJr+6rRhU40somIyA+lbrIluQeeCUD1s7fTtvSLqGeI5jDz+UCJMWY68BdgJ+B0YAru9aZLjDGvAjXALOCZcG9KREQ8kLXtfnRUL6Vh7jNU/utaSk67hkB2ftQ+P5qj+OqBg4EjgVrc03tnWGtfwz2l9zDwGlCB23M6LVrZRERk7Ybtewqpo7amp6WB5Y9cQ09Ha9Q+O6o36lprPwR2X8v2TmBG+CEiIjHC8fnJnzKd8tm/o6NqCVVP3kTBUefhOAPfv9FURyIisl7+1AwKjv4dvmAaLfY9Vr7xaFQ+VwVKRER+VHJuKflTprtz9r3xCM0L3x3wz1SBEhGRiKSN2Zah+5wEQNWTN9FR9e2Afp4KlIiIRCx7p8PJ2HJPQp1tVD5yNd0tAzfYWgVKREQi5jgOuQf/jGDRGLrqq6h5/q4B+ywVKBER6RNfIEjB1AtJzh+JPy17wD4nqsPMRUQkMSRlDaP0zOsH9DPUgxIRkZikAiUiIjFJBUpERGKSCpSIiMQkFSgREYlJKlAiIhKTVKBERCQmJcJ9UH6AyspKr3OIiEgf9fru9q+5LxEKVBHAiSee6HUOERHZcEXA1703JEKBmgvsgbsSb7fHWUREpG/8uMVp7po7nFAoFP04IiIiP0KDJEREJCapQImISExSgRIRkZikAiUiIjFJBUpERGKSCpSIiMQkFSgREYlJKlAiIhKTEmEmiY1ijNkauA2YAHwDnGat/Z87muOFMWY/4GpgU6AKuM5ae7sxJhn4GzAVd8aNG6y1V3mXdMMZY3KAj4E/WGtnJ1jbioBbgZ8AbcAd1toZidBGY8zOwF8BA1QDV1tr74r3thljdgSettbmh5+vtz3GmGOAK3FnT3gNmGatrYp68AitpX35wF+ASYADPAeca61dGd7fb+0b1D2o8D+k/wAPAznALOAFY0yWp8E2kDFmOPBv4Arc9hwPXGWMOQC4DPeLYQywA3CKMeZkr7JupNuAkl7PE6lt/8GdtqsA2Bm3LScQ5200xvhw2/ZXa2027r/Nv4X/QIzLthljHGPMGcALQHKvXetsjzFmC+BuYBowDPgKeCiKsSO2nvbdBXQBo3D/EB4C3Bw+pl/bN6gLFLA3ELDW3mit7bTWPgR8BhzrbawNNhJ4wFr7uLW2J9wTfBXYDTgFmGWtXWmtXQL8CTjbq6AbyhhzCpAFfNJrc6K0bSdgNPB/1to2a+1i3H+jrxD/bRwC5AOOMcYBQrhfch3Eb9suA36O+wdhb+trz0nAU9baN621bcBFwG7GmE2jlLkv/qd94T80eoDLrLXN1to64E5g9/BL+rV9g71AbQF8sca2hcBWHmTZaNbaN6y1P1v13BgzFHci3Y9wu9uf93p53LXTGDMKmAmc1mtbDgnQtrCJuIX3UmPMMmPM18AUoJU4b6O1dgXuaa97gE7ciUEvxu0txmvbbrPWTgQ+WLUhgn+PW/TeZ61tAZYSm+39n/aF//CdbK1d1Ot1k3G/Y6Cf2zfYr0FlAC1rbGsB0jzI0q+MMdnAk8B7wLzw5t5tjat2GmP8wH3AedbaSmPMql0Z4Z9x27ZeVv1B8RpuT2ocMAf3eg3EcRvDf3m3ASfgnobeFXgMqAu/JO7aZq0tX8vmH/v3GDffOeto3w8YY87DLVC7hjf1a/sGe4FqBlLX2JYGNHmQpd8YYzbDPd//OXAi37exd1vjrZ0zAGutfWyN7c3hn/HctlXagQZr7aXh5wuMMXfhnjKC+G7jkcBu1trzw89fM8bcTWK0rbcf+/eYEN85xpgAcBNwGLCPtXZheFe/tm+wn+L7HPdiZm/j+GH3PK4YY/bE7TU9AUwNX8tYCVTyw7bGWzuPA6YaY+qMMXW4pwxuwR3YEu9tW2UhkBYevLNKEpAI//2GA8E1tnXh9g7jvW2rRfD/2g++c4wxacAI4qi9xphM4EXcASA7Wmvn99rdr+0b7D2oV3Av2v4G9/z4UbjDzR/3NNUGMsaMAZ4GLrHW3rTG7nuBmcaYj3G74efhDhWNC9bacb2fG2PmAzeGh5k3Ecdt6+VF3C/s640xv8X9H/103AvV3xDfbXwBd0TpWbgX1bcDzgTOAL4jvtu2pvX9v/YA8KYxZm/gHeAq4CNr7ZdeBN1AD+F2bvYIX2PqrV/bN6h7UNbaDuAg3MJUC1wCTLbWVq/3wNh1DpCJ+0XQ1OtxDfAH4FPcUYpzca8D3OZd1H6VEG0Lj3raC/f6UwXu9adrrbX/Js7baK39DPc039m4150eAH5nrf0Pcd62tVhne6y1n+AO8rkNqAHGA0d7E7PvjDETgIOBHYGqXt8xZdD/7dOKuiIiEpMGdQ9KRERilwqUiIjEJBUoERGJSSpQIiISk1SgREQkJqlAiYhITBrsN+qKbDBjzGy+n6pnbS7DnU3+FSDTWhuV6WzC8xa+BZy8vhskw/PjvQv81Fpro5FNpC/UgxLZcOfizlxdhLssBrg3MK7a9ifg7fDvzWs5fqD8H7Dgx+7et9b2AJcT3zfFSgLTjboi/cAYsyXuUhmjwmsAeZUjBXfqoH2stZ9GeMzXwOnW2lcHMptIX+kUn8gACs9JtvoUnzEmhLua7EW4c+19gLvI2/nAT4EG4CJr7b3h4zOB63GXDw8BL+Mur72upRCOA+p6FydjzAzgLCAPd/2zi621z/U65nHc3uCr/dBkkX6jU3wi0Xc18GvcJd1HAB/iFqYdcNdIut0Ys2pdoTtwC9kBuPP0hYDnjTHr+uPyENw5/AAwxkwJf9ZJuLNqPwM8aozJ6nXMHGDf9byniCdUoESi72Zr7SvhZQqexl0r5+LwQIUbcNfTGWWMGY3bIzrBWjs33Cv6KTASOHAd77097iSlq4zEXWfq2/Cpx8txJ23t7PWaz3Fn3f7BjPEiXtNfTCLR13u57BZgibV21cXgtvDPILBJ+HfbawVhcBeAM7jFbU0FuLNIr3If7kjDb4wx83BXWf6Htba112tWhH/m97EdIgNKPSiR6Otc43nPOl6XFH7ttsA2vR6bAf9YxzE9gLPqSXjpmIm4Pa63gWnAx+FBHaus+h7ojrgFIlGgAiUSu74AAkC6tXaRtXYR7jpR1+EWqbWpxB0MAYAx5kjgbGvtC9bac3F7Xo24a/qsktfrWJGYoVN8IjHKWmuNMU8C/zTGnIO72u4s3MEVC9dx2Dxg617P/cB1xpjluCMGdwYKw7+vsjXusvK9Tz2KeE49KJHYdgpuMXkCd3XWbGA/a23dOl7/DO5oPwCstY8CM3F7XV8CVwC/tNa+3OuYPYE51lqd4pOYoht1RRKIMSYNWAIcaK39MILX+4BvcUcKvjHA8UT6RD0okQRirW3B7S2dE+EhRwDfqDhJLFKBEkk8fwYmmDXGpq8p3Hu6BPhZVFKJ9JFO8YmISExSD0pERGKSCpSIiMQkFSgREYlJKlAiIhKTVKBERCQm/T8DXY63ypXELAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_y(results):\n",
    "    plot(results.y, color='C1', label='y')\n",
    "\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Length (m)')\n",
    "    \n",
    "plot_y(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
